Fit length models

# Load libraries
library(tidyverse)
library(tidylog)
library(sdmTMB)
library(patchwork)
library(viridis)
library(RColorBrewer)
library(modelr)
library(ggstats)
library(ggspatial)
library(ggsidekick)
theme_set(theme_sleek())

home <- here::here()
source(paste0(home, "/R/functions/map-plot.R"))

remove_attributes <- function(x) {attributes(x) <- NULL; return(x)}

Read data and prediction grid, scale variables

# Read data
d <- readr::read_csv(paste0(home, "/data/clean/larval_size.csv")) |>
  drop_na(temp) |>
  drop_na(chl) |>
  mutate(
    yday_ct = yday - mean(yday),
    year_f = as.factor(year),
    species_f = as.factor(species),
    year_ct = year - median(year)
  ) |>
  mutate(across(
    .cols = c("temp", "chl"),
    .fns = scale,
    .names = "{.col}_sc"
  )) |> 
  mutate(temp_sq = temp_sc*temp_sc) |> 
  filter(!species == "Anguilla anguilla") |> 
  filter(!species == "Sardina pilchardus") |> 
  as.data.frame() |> 
  mutate(across(everything(), as.vector))

#d <- lapply(d, remove_attributes)

coul <- brewer.pal(11, "Spectral")
coul <- colorRampPalette(coul)(length(unique(d$species)))

ggplot(d, aes(length_mm, fill = species)) +
  geom_histogram() +
  scale_fill_manual(values = coul, name = "Species") +
  coord_cartesian(expand = 0) +
  labs(y = "Count", x = "Length (mm)", tag = "b)") +
  theme(
    legend.text = element_text(face = "italic", size = 7),
    legend.key.size = unit(0.25, "cm"),
    legend.position.inside = c(0.2, 0.69),
    plot.tag = element_text()
  )

# Load prediction grid
pred_grid <- readr::read_csv(paste0(home, "/data/clean/pred_grid.csv")) |>
  drop_na(temp) |>
  drop_na(chl) |>
  filter(year %in% unique(d$year)) |>
  mutate(
    temp_sc = (temp - mean(d$temp, na.rm = TRUE)) / sd(d$temp, na.rm = TRUE),
    temp_sq = temp_sc*temp_sc,
    chl_sc = (chl - mean(chl, na.rm = TRUE)) / sd(chl, na.rm = TRUE),
    depth_sc = as.numeric(scale(depth)),
    year_f = as.factor(year),
    year_ct = 0,
    yday_ct = 0
  ) |>
  mutate(keep = ifelse(lon < 10 & lat < 57.15, "N", "Y")) |>
  filter(keep == "Y") |>
  dplyr::select(-keep)

# Read density data
b <- readr::read_csv(paste0(home, "/data/clean/larval_density.csv")) |>
  drop_na(no_m2) |> 
  mutate(year_f  = as.factor(year),
         depth_sc = (depth - mean(pred_grid$depth, na.rm = TRUE))/sd(pred_grid$depth, na.rm = TRUE)) |> 
  mutate(sum_no_m2 = sum(no_m2), .by = c(year, species)) |> 
  filter(sum_no_m2 > 0)

length(unique(b$haul_id))

[1] 1530

Plot covariates

# First gather the elements common to all the covariate plots
covMapTheme <- plot_map +
  theme(
    legend.position.inside = c(0.25, 0.11),
    legend.direction = "horizontal",
    legend.key.width = unit(0.35, "cm"), # Decreased to fit multicovariates
    legend.key.height = unit(0.2, "cm"),
    legend.text = element_text(size = 7),
    legend.title = element_text(size = 8),
    plot.tag = element_text()
  ) +
  geom_sf() +
  annotate("text",
    label = "Sweden", x = xmin2 + 0.95 * xrange, y = ymin2 + 0.75 * yrange,
    color = "gray50", size = 2.5
  ) +
  annotate("text",
    label = "Norway", x = xmin2 + 0.08 * xrange, y = ymin2 + 0.95 * yrange,
    color = "gray50", size = 2.5
  ) +
  annotate("text",
    label = "Denmark", x = xmin2 + 0.42 * xrange, y = ymin2 + 0.45 * yrange,
    color = "gray50", size = 2.5
  ) +
  guides(fill = guide_colorbar(
    position = "inside",
    title.position = "top",
    title.hjust = 0.5
  )) +
  NULL

# Generate the data and covariate maps
nb.cols <- length(unique(d$species))
mycolors <- colorRampPalette(brewer.pal(12, "Set3"))(nb.cols)

p_dat <- covMapTheme +
  geom_point(
    data = d, aes(X * 1000, Y * 1000, color = species),
    size = 0.5, alpha = 1
  ) +
  scale_color_manual(values = mycolors) +
  labs(tag = "a)") +
  theme(axis.title.x = element_blank()) +
  annotation_scale(bar_cols = c("grey30", "white"), height = unit(0.1, "cm")) +
  annotation_north_arrow(
    location = "tr", which_north = "true", height = unit(0.85, "cm"),
    width = unit(0.85, "cm"), pad_x = unit(0.1, "in"), pad_y = unit(0.1, "in"),
    style = north_arrow_fancy_orienteering(
      fill = c("grey40", "white"),
      line_col = "grey20"
    )
  )

p_tempCovMap <- covMapTheme +
  geom_raster(
    data = pred_grid |>
      group_by(X, Y) |>
      summarise(mean_temp = mean(temp)),
    aes(X * 1000, Y * 1000, fill = mean_temp)
  ) +
  labs(tag = "b)") +
  scale_fill_viridis(option = "magma", name = "Temperature (°C)") +
  theme(
    axis.text.y = element_blank(),
    axis.title.y = element_blank()
  ) +
  geom_sf()

group_by: 2 grouping variables (X, Y)

summarise: now 4,689 rows and 3 columns, one group variable remaining (X)

skag_kat_border <- data.frame(lat1 = 57.75420, 
                              lat2 = 57.91694673025835,
                              lon1 = 10.61630,
                              lon2 = 11.502651981474733) 

skag_kat_border <- skag_kat_border |> 
  add_utm_columns(utm_names = c("X1", "Y1"),
                  ll_names = c("lon1", "lat1")) |> 
  add_utm_columns(utm_names = c("X2", "Y2"),
                  ll_names = c("lon2", "lat2"))

Detected UTM zone 32N; CRS = 32632.

Visit https://epsg.io/32632 to verify.

Detected UTM zone 32N; CRS = 32632.

Visit https://epsg.io/32632 to verify.

ann_text <- data.frame(year = 2012, 
                       label = c("Skagerrak", "Kattegat"),
                       lon = c(10.55, 11.2),
                       lat = c(58.1, 57.5)) |> 
  add_utm_columns(ll_names = c("lon", "lat"))

Detected UTM zone 32N; CRS = 32632.
Visit https://epsg.io/32632 to verify.

p_tempCovMap <- p_tempCovMap + 
  geom_segment(data = skag_kat_border, aes(x = X1*1000, xend = X2*1000,
                                           y = Y1*1000, yend = Y2*1000),
               alpha = 0.75, linetype = 3, color = "white") + 
  geom_text(data = ann_text, aes(x = X*1000, y = Y*1000, label = label),
            color = "white", size = 2.7)
  
  
p_chlCovMap <- covMapTheme +
  geom_raster(
    data = pred_grid |>
      group_by(X, Y) |>
      summarise(mean_chl = mean(chl)),
    aes(X * 1000, Y * 1000, fill = mean_chl)
  ) +
  labs(tag = "c)") +
  scale_fill_viridis(name = expression("Chl (mg/l"^3 * ")")) +
  theme(
    axis.text.y = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  ) +
  geom_sf()

group_by: 2 grouping variables (X, Y)

summarise: now 4,689 rows and 3 columns, one group variable remaining (X)

# Plot over time
trends <- pred_grid |>
  pivot_longer(c(temp, chl), names_to = "var") |>
  mutate(var = fct_recode(var,
    "Temperature" = "temp",
    "Chlorophyll" = "chl"
  )) |>
  summarise(value = mean(value), .by = c(var, year))

pivot_longer: reorganized (temp, chl) into (var, value) [was 135981x15, now 271962x15]

mutate: converted 'var' from character to factor (0 new NA)

summarise: now 58 rows and 3 columns, ungrouped

ct <- trends |>
  filter(var == "Chlorophyll") |>
  ggplot(aes(year, value)) +
  geom_point(color = "grey30") +
  geom_smooth(method = "lm", color = "steelblue") +
  labs(x = "Year", y = expression("Chl (mg/l"^3 * ")"), tag = "d)")

filter: removed 29 rows (50%), 29 rows remaining

tt <- trends |>
  filter(var == "Temperature") |>
  ggplot(aes(year, value)) +
  geom_point(color = "grey30") +
  geom_smooth(method = "lm", color = "steelblue") +
  labs(x = "Year", y = "Temperature (°C)", tag = "e)")

filter: removed 29 rows (50%), 29 rows remaining

# Combined plot
(p_dat + p_tempCovMap + p_chlCovMap) / free(tt + ct) +
  plot_layout(heights = c(0.9, 1))

Warning: `guide_colourbar()` needs continuous scales.

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

# Save
ggsave(paste0(home, "/figures/data_map.pdf"), width = 21, height = 18, units = "cm")

Warning: `guide_colourbar()` needs continuous scales.

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Fit models by species

  # db |>
  #   mutate(pres = ifelse(no_m2 > 0, 1, 0)) |>
  #   ggplot(aes(X, Y, color = pres)) +
  #   geom_point() +
  #   facet_wrap(~year)
  # hist(db$no_m2)
  # summary(db$no_m2)
  # db |>
  #   ggplot(aes(X, Y, color = no_m2)) +
  #   scale_color_viridis(trans = "sqrt") +
  #   geom_point() +
  #   facet_wrap(~year)
  # dd |>
  #   ggplot(aes(X, Y)) +
  #   geom_point() +
  #   facet_wrap(~year)

pars_list <- list()
res_list <- list()
res_dens_list <- list()
spatial_list <- list()
spatial_density_list <- list()
range_list <- list()
genq_list <- list()
index_list <- list()
dens_index_list <- list()

for (i in unique(d$species)) {
  print(i)

  dd <- d |> filter(species == i)
  db <- b |> filter(species == i)

  if (unique(dd$species) %in% c(
    "Ammodytidae", "Crystallogobius linearis",
    "Syngnathus rostellatus", "Aphia minuta"
  )) {

    mesh <- make_mesh(dd, xy_cols = c("X", "Y"), cutoff = 8)
    mesh_dens <- make_mesh(db, xy_cols = c("X", "Y"), cutoff = 8)

    ggplot() +
      inlabru::gg(mesh$mesh) +
      coord_fixed() +
      geom_point(aes(X, Y), data = dd, alpha = 0.2, size = 0.5) +
      annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.1, vjust = 2) +
      labs(x = "Easting (km)", y = "Northing (km)")

    ggsave(paste0(home, paste("/figures/supp/size_mesh_", i, "_.pdf")), width = 14, height = 14, units = "cm")

    mlength <- sdmTMB(length_mm ~ temp_sc + chl_sc + year_f + yday_ct,
      data = dd,
      mesh = mesh,
      family = gengamma(link = "log"),
      spatiotemporal = "off",
      spatial = "on",
      time = "year")

    mdens <- sdmTMB(no_m2 ~ year_f,
      data = db,
      mesh = mesh_dens,
      family = delta_gamma(type = "poisson-link"),
      spatiotemporal = "iid",
      spatial = "on",
      time = "year")
    
    mlength <- run_extra_optimization(mlength)
  } else {

    mesh <- make_mesh(dd, xy_cols = c("X", "Y"), cutoff = 4)
    mesh_dens <- make_mesh(db, xy_cols = c("X", "Y"), cutoff = 4)
  
    ggplot() +
      inlabru::gg(mesh$mesh) +
      coord_fixed() +
      geom_point(aes(X, Y), data = dd, alpha = 0.2, size = 0.5) +
      annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.1, vjust = 2) +
      labs(x = "Easting (km)", y = "Northing (km)")

    ggsave(paste0(home, paste("/figures/supp/size_mesh_", i, "_.pdf")), width = 14, height = 14, units = "cm")

    mlength <- sdmTMB(length_mm ~ temp_sc + chl_sc + year_f + yday_ct,
      data = dd,
      mesh = mesh,
      family = gengamma(link = "log"),
      spatiotemporal = "off",
      spatial = "on",
      time = "year")
    
    mlength <- run_extra_optimization(mlength)
    
    mdens <- sdmTMB(no_m2 ~ year_f,
      data = db,
      mesh = mesh_dens,
      family = delta_gamma(type = "poisson-link"),
      spatiotemporal = "iid",
      spatial = "on",
      time = "year")

  }

  print(mlength)
  sanity(mlength)
  
  print(mdens)
  sanity(mdens)
  
  # Residuals
  # For the density model we'll need to use dharma residuals because 2 models (hurdle), so lets do that for length also
  res <- simulate(mlength, nsim = 500, type = "mle-mvn") |>
    dharma_residuals(mlength)
  dd$observed <- res$observed
  dd$expected <- res$expected
  
  res_dens <- simulate(mdens, nsim = 500, type = "mle-mvn") |>
    dharma_residuals(mdens)
  db$observed <- res_dens$observed
  db$expected <- res_dens$expected
  
  res_list[[i]] <- dd
  res_dens_list[[i]] <- db
  
  # Get spatial predictions
  p <- predict(mlength, newdata = pred_grid |>
                 filter(year %in% unique(dd$year_f)))

  # Predict with density model
  pdens <- predict(mdens, 
                   newdata = pred_grid |>
                     filter(year %in% unique(db$year_f))) |> 
    mutate(pred_dens = exp(est1 + est2)) |> 
    dplyr::select(-est_non_rf1, -est_non_rf2, -est_rf2, -est_rf2,
                  -omega_s1, -omega_s2, #-epsilon_st1, -epsilon_st2,
                  -est1, -est2)
  
  spatial_list[[i]] <- p |> mutate(species = i)
  spatial_density_list[[i]] <- pdens |> mutate(species = i)
  
  # Get index - length
  p_ind <- predict(mlength,
    newdata = pdens |> filter(year %in% unique(dd$year)),
    return_tmb_object = TRUE
  )

  ind <- get_index(p_ind,
    area = 1 / nrow(pred_grid |> filter(year == max(year))),
    bias_correct = TRUE
  ) |> 
    mutate(type = "unweighted")
  
  # Get index - density
  p_dens <- predict(mdens, 
                    newdata = pred_grid |>
                      filter(year %in% unique(db$year_f)),
                    return_tmb_object = TRUE) 
  
  d_ind <- get_index(p_dens, area = 9*1e6, bias_correct = TRUE)
  
  dens_index_list[[i]] <- d_ind |> mutate(species = i)
  
  # Weighted index
  weight_sum <- pdens |> 
    summarise(pred_dens_sum = sum(pred_dens), .by = year)
  
  ind_w <- get_index(p_ind,
    area = p_ind$data$pred_dens,
    bias_correct = TRUE
  ) |> 
    left_join(weight_sum) |> 
    mutate(est = est / pred_dens_sum,
           lwr = lwr / pred_dens_sum,
           upr = upr / pred_dens_sum)
  
  ind <- ind |> bind_rows(ind_w |> mutate(type = "weighted"))
  
  index_list[[i]] <- ind |> mutate(species = i)
  
  # Get ranges
  range_length <- tidy(mlength, effects = "ran_pars") |>
    filter(term == "range") |> 
    mutate(cutoff = ifelse(i %in% c(
      "Ammodytidae", "Crystallogobius linearis",
      "Syngnathus rostellatus", "Aphia minuta"
    ),
    8, 4
    )) |> 
    mutate(model = "length")
  
  range_dens <- tidy(mdens, effects = "ran_pars") |>
    filter(term == "range") |> 
    mutate(cutoff = ifelse(i %in% c(
      "Ammodytidae", "Crystallogobius linearis",
      "Syngnathus rostellatus", "Aphia minuta"
    ),
    8, 4
    )) |> 
    mutate(model = "density")

  range <- bind_rows(range_length, range_dens)
  
  # Get gengamma Q and phi
  genq_list[[i]] <- tibble(
    term = "gengamma_Q",
    estimate = mlength$sd_report$par.fixed["gengamma_Q"]
  ) |>
    bind_rows(tidy(mlength, effects = "ran_pars") |> filter(term == "phi")) |>
    mutate(species = i)

  range_list[[i]] <- range |> mutate(species = i)
  
  # Other parameters
  pars_list[[i]] <- tidy(mlength) |> mutate(species = i)
  
}

[1] "Aphia minuta"

filter: removed 37,448 rows (86%), 6,325 rows remaining

filter: removed 16,734 rows (92%), 1,530 rows remaining

Loading required namespace: INLA

Spatial model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc + chl_sc + year_f + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: gengamma(link = 'log')
 
            coef.est coef.se
(Intercept)     3.45    0.02
temp_sc        -0.07    0.01
chl_sc         -0.01    0.01
year_f1994     -0.12    0.03
year_f1995      0.09    0.04
year_f1996     -0.13    0.05
year_f1997     -0.14    0.04
year_f1998      0.14    0.02
year_f1999     -0.06    0.02
year_f2000      0.11    0.02
year_f2001      0.17    0.03
year_f2002      0.04    0.02
year_f2003     -0.13    0.03
year_f2004      0.13    0.02
year_f2005      0.05    0.02
year_f2006     -0.07    0.02
year_f2007      0.20    0.03
year_f2008      0.17    0.03
year_f2009      0.01    0.03
year_f2010     -0.12    0.03
year_f2012      0.10    0.03
year_f2013     -0.03    0.03
year_f2014      0.16    0.03
year_f2015      0.09    0.03
year_f2016      0.03    0.03
year_f2017     -0.01    0.04
year_f2018      0.01    0.03
year_f2019      0.08    0.04
year_f2020      0.05    0.04
year_f2021     -0.04    0.03
year_f2022      0.09    0.03
yday_ct         0.00    0.00

Dispersion parameter: 0.17
Generalized gamma Q: 1.14
Matérn range: 19.62
Spatial SD: 0.08
ML criterion at convergence: 21545.082

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: no_m2 ~ year_f
Mesh: mesh_dens (isotropic covariance)
Time column: year
Data: db
Family: delta_gamma(link1 = 'log', link2 = 'log', type = 'poisson-link')

Delta/hurdle model 1: -----------------------------------
Family: binomial(link = 'log') 
            coef.est coef.se
(Intercept)    -0.35    0.43
year_f1994     -0.51    0.52
year_f1995      0.12    0.50
year_f1996     -1.51    0.58
year_f1997     -0.99    0.55
year_f1998      0.42    0.50
year_f1999      1.18    0.49
year_f2000      0.47    0.48
year_f2001      0.87    0.47
year_f2002     -0.43    0.50
year_f2003     -0.18    0.50
year_f2004     -0.02    0.48
year_f2005      1.33    0.52
year_f2006      0.58    0.48
year_f2007      1.87    0.50
year_f2008      0.27    0.49
year_f2009      0.23    0.48
year_f2010     -0.29    0.50
year_f2012      1.13    0.49
year_f2013      1.24    0.50
year_f2014      1.51    0.50
year_f2015      0.27    0.47
year_f2016      0.31    0.47
year_f2017     -0.63    0.49
year_f2018      0.61    0.48
year_f2019     -1.11    0.52
year_f2020     -0.25    0.52
year_f2021      0.15    0.50
year_f2022     -0.73    0.54

Matérn range: 77.22
Spatial SD: 0.66
Spatiotemporal IID SD: 0.67

Delta/hurdle model 2: -----------------------------------
Family: Gamma(link = 'log') 
            coef.est coef.se
(Intercept)    -4.56    0.26
year_f1994      0.09    0.35
year_f1995      1.41    0.36
year_f1996     -0.06    0.39
year_f1997      0.51    0.37
year_f1998      0.36    0.35
year_f1999      0.19    0.35
year_f2000      0.40    0.32
year_f2001      0.27    0.31
year_f2002      0.86    0.33
year_f2003      0.12    0.32
year_f2004      0.56    0.33
year_f2005      0.32    0.36
year_f2006      0.09    0.33
year_f2007      0.10    0.36
year_f2008      0.05    0.33
year_f2009      0.12    0.32
year_f2010      0.22    0.33
year_f2012     -0.54    0.34
year_f2013      0.28    0.35
year_f2014     -0.12    0.36
year_f2015     -0.11    0.31
year_f2016      0.16    0.31
year_f2017     -0.21    0.32
year_f2018      0.28    0.32
year_f2019      0.03    0.34
year_f2020     -0.36    0.34
year_f2021     -0.21    0.34
year_f2022      0.41    0.39

Dispersion parameter: 3.36
Matérn range: 19.98
Spatial SD: 0.34
Spatiotemporal IID SD: 0.65

ML criterion at convergence: -1288.014

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameters don't look unreasonably large

Simulating ■■■■                               9% | ETA: 10s

Simulating ■■■■■■                            17% | ETA:  9s

Simulating ■■■■■■■■■■■■■■                    44% | ETA:  6s

Simulating ■■■■■■■■■■■■■■■■■■■■■■            71% | ETA:  3s

Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■    98% | ETA:  0s

Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

filter: no rows removed

filter: no rows removed

mutate: new variable 'pred_dens' (double) with 135,981 unique values and 0% NA

mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA

filter: no rows removed

filter: removed 131,292 rows (97%), 4,689 rows remaining

mutate: changed 29 values (100%) of 'type' (0 new NA)

filter: no rows removed

mutate: new variable 'species' (character) with one unique value and 0% NA

summarise: now 29 rows and 2 columns, ungrouped

Joining with `by = join_by(year)`
left_join: added one column (pred_dens_sum)
> rows only in x 0
> rows only in y ( 0)
> matched rows 29
> ====
> rows total 29
mutate: changed 29 values (100%) of 'est' (0 new NA)
changed 29 values (100%) of 'lwr' (0 new NA)
changed 29 values (100%) of 'upr' (0 new NA)
mutate: changed 29 values (100%) of 'type' (0 new NA)
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA

[1] "Clupea harengus"

filter: removed 32,194 rows (74%), 11,579 rows remaining
filter: removed 16,734 rows (92%), 1,530 rows remaining

Spatial model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc + chl_sc + year_f + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: gengamma(link = 'log')
 
            coef.est coef.se
(Intercept)     3.34    0.01
temp_sc        -0.05    0.01
chl_sc          0.02    0.01
year_f1994     -0.05    0.01
year_f1995     -0.05    0.01
year_f1996     -0.06    0.03
year_f1997     -0.12    0.02
year_f1998      0.03    0.02
year_f1999     -0.01    0.01
year_f2000      0.01    0.01
year_f2001      0.17    0.02
year_f2002      0.00    0.01
year_f2003     -0.13    0.03
year_f2004      0.07    0.02
year_f2005      0.04    0.02
year_f2006     -0.02    0.02
year_f2007      0.02    0.02
year_f2008      0.01    0.02
year_f2009     -0.03    0.01
year_f2012     -0.23    0.02
year_f2013     -0.02    0.02
year_f2014     -0.05    0.02
year_f2015     -0.08    0.02
year_f2016      0.02    0.03
year_f2017     -0.11    0.02
year_f2018     -0.09    0.02
year_f2019     -0.20    0.01
year_f2020     -0.07    0.02
year_f2021     -0.08    0.02
year_f2022     -0.15    0.02
yday_ct         0.00    0.00

Dispersion parameter: 0.12
Generalized gamma Q: 0.56
Matérn range: 10.40
Spatial SD: 0.08
ML criterion at convergence: 30631.810

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: no_m2 ~ year_f
Mesh: mesh_dens (isotropic covariance)
Time column: year
Data: db
Family: delta_gamma(link1 = 'log', link2 = 'log', type = 'poisson-link')

Delta/hurdle model 1: -----------------------------------
Family: binomial(link = 'log') 
            coef.est coef.se
(Intercept)     4.60  355.94
year_f1994     -2.95  355.94
year_f1995     -2.54  355.94
year_f1996     -6.37  355.94
year_f1997     -4.42  355.94
year_f1998     -4.77  355.94
year_f1999     -3.15  355.94
year_f2000     -3.93  355.94
year_f2001     -4.56  355.94
year_f2002     -3.19  355.94
year_f2003     -6.32  355.94
year_f2004     -5.18  355.94
year_f2005     -4.45  355.94
year_f2006     -5.14  355.94
year_f2007     -4.14  355.94
year_f2008     -4.89  355.94
year_f2009     -4.04  355.94
year_f2010     -7.77  355.94
year_f2012     -3.44  355.94
year_f2013     -5.25  355.94
year_f2014     -4.20  355.94
year_f2015     -4.80  355.94
year_f2016     -5.78  355.94
year_f2017     -4.12  355.94
year_f2018     -5.06  355.94
year_f2019     -4.23  355.94
year_f2020     -3.51  355.94
year_f2021     -5.75  355.94
year_f2022     -3.54  355.94

Matérn range: 66.87
Spatial SD: 0.58
Spatiotemporal IID SD: 0.68

Delta/hurdle model 2: -----------------------------------
Family: Gamma(link = 'log') 
            coef.est coef.se
(Intercept)    -6.42  355.94
year_f1994      2.60  355.94
year_f1995      2.34  355.94
year_f1996      2.96  355.94
year_f1997      2.60  355.94
year_f1998      2.66  355.94
year_f1999      2.69  355.94
year_f2000      2.54  355.94
year_f2001      2.69  355.94
year_f2002      2.42  355.94
year_f2003      1.92  355.94
year_f2004      2.23  355.94
year_f2005      2.18  355.94
year_f2006      1.86  355.94
year_f2007      2.26  355.94
year_f2008      1.81  355.94
year_f2009      2.01  355.94
year_f2010      1.75  355.94
year_f2012      1.65  355.94
year_f2013      1.58  355.94
year_f2014      1.63  355.94
year_f2015      1.44  355.94
year_f2016      1.92  355.94
year_f2017      1.68  355.94
year_f2018      2.10  355.94
year_f2019      1.84  355.94
year_f2020      1.15  355.94
year_f2021      2.25  355.94
year_f2022      1.49  355.94

Dispersion parameter: 2.62
Matérn range: 65.72
Spatial SD: 0.54
Spatiotemporal IID SD: 0.26

ML criterion at convergence: -1215.562

See ?tidy.sdmTMB to extract these values as a data frame.

**Possible issues detected! Check output of sanity().**

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✖ `b_j` standard error may be large
ℹ Try simplifying the model, adjusting the mesh, or adding priors

✖ `b_j2` standard error may be large
ℹ Try simplifying the model, adjusting the mesh, or adding priors

✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameters don't look unreasonably large
Simulating ■■■■■■■■■■■■                      35% | ETA:  2s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

Simulating ■■                                 3% | ETA: 31s
Simulating ■■■                                6% | ETA: 31s
Simulating ■■■■■                             15% | ETA: 28s
Simulating ■■■■■■■■                          24% | ETA: 26s
Simulating ■■■■■■■■■■■                       33% | ETA: 23s
Simulating ■■■■■■■■■■■■■                     42% | ETA: 20s
Simulating ■■■■■■■■■■■■■■■■                  50% | ETA: 17s
Simulating ■■■■■■■■■■■■■■■■■■■               59% | ETA: 14s
Simulating ■■■■■■■■■■■■■■■■■■■■■             68% | ETA: 11s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■          77% | ETA:  8s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■       86% | ETA:  5s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■     95% | ETA:  2s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

filter: removed 4,689 rows (3%), 131,292 rows remaining
filter: no rows removed
mutate: new variable 'pred_dens' (double) with 135,981 unique values and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 4,689 rows (3%), 131,292 rows remaining
filter: removed 131,292 rows (97%), 4,689 rows remaining
mutate: changed 28 values (100%) of 'type' (0 new NA)
filter: no rows removed
mutate: new variable 'species' (character) with one unique value and 0% NA
summarise: now 29 rows and 2 columns, ungrouped
Joining with `by = join_by(year)`left_join: added one column (pred_dens_sum)
           > rows only in x    0
           > rows only in y  ( 1)
           > matched rows     28
           >                 ====
           > rows total       28
mutate: changed 28 values (100%) of 'est' (0 new NA)
        changed 28 values (100%) of 'lwr' (0 new NA)
        changed 28 values (100%) of 'upr' (0 new NA)
mutate: changed 28 values (100%) of 'type' (0 new NA)
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA

[1] "Crystallogobius linearis"

filter: removed 29,637 rows (68%), 14,136 rows remaining
filter: removed 16,734 rows (92%), 1,530 rows remaining

Spatial model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc + chl_sc + year_f + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: gengamma(link = 'log')
 
            coef.est coef.se
(Intercept)     3.09    0.03
temp_sc         0.00    0.01
chl_sc          0.00    0.01
year_f1994      0.06    0.02
year_f1995      0.19    0.02
year_f1996      0.15    0.03
year_f1997      0.15    0.03
year_f1998      0.16    0.02
year_f1999      0.05    0.02
year_f2000      0.13    0.02
year_f2001      0.10    0.03
year_f2002      0.08    0.03
year_f2003      0.03    0.03
year_f2004      0.04    0.02
year_f2005      0.00    0.02
year_f2006      0.06    0.02
year_f2007      0.00    0.03
year_f2008      0.02    0.02
year_f2009      0.05    0.02
year_f2010      0.04    0.04
year_f2012      0.00    0.03
year_f2013     -0.02    0.03
year_f2014      0.11    0.03
year_f2015     -0.10    0.03
year_f2016      0.02    0.03
year_f2017     -0.12    0.03
year_f2018     -0.02    0.03
year_f2019     -0.06    0.02
year_f2020     -0.18    0.03
year_f2021     -0.13    0.02
year_f2022     -0.04    0.03
yday_ct         0.00    0.00

Dispersion parameter: 0.23
Generalized gamma Q: 0.04
Matérn range: 61.94
Spatial SD: 0.09
ML criterion at convergence: 42417.797

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: no_m2 ~ year_f
Mesh: mesh_dens (isotropic covariance)
Time column: year
Data: db
Family: delta_gamma(link1 = 'log', link2 = 'log', type = 'poisson-link')

Delta/hurdle model 1: -----------------------------------
Family: binomial(link = 'log') 
            coef.est coef.se
(Intercept)     0.44    0.62
year_f1994      0.01    0.42
year_f1995      0.35    0.42
year_f1996     -0.76    0.42
year_f1997     -0.40    0.41
year_f1998      0.19    0.40
year_f1999      0.56    0.41
year_f2000     -0.05    0.41
year_f2001      0.22    0.39
year_f2002     -0.11    0.41
year_f2003     -0.53    0.40
year_f2004      0.08    0.40
year_f2005      0.52    0.44
year_f2006     -0.26    0.41
year_f2007     -0.04    0.39
year_f2008      0.31    0.41
year_f2009      0.09    0.40
year_f2010     -1.18    0.41
year_f2012      0.27    0.40
year_f2013      0.14    0.40
year_f2014     -0.04    0.40
year_f2015     -0.41    0.39
year_f2016      0.26    0.41
year_f2017      0.52    0.42
year_f2018      0.73    0.49
year_f2019     -0.16    0.40
year_f2020     -0.32    0.42
year_f2021      0.30    0.45
year_f2022     -0.26    0.44

Matérn range: 129.39
Spatial SD: 0.96
Spatiotemporal IID SD: 0.31

Delta/hurdle model 2: -----------------------------------
Family: Gamma(link = 'log') 
            coef.est coef.se
(Intercept)    -4.64    0.28
year_f1994      0.84    0.36
year_f1995      0.83    0.36
year_f1996      0.22    0.35
year_f1997      0.42    0.34
year_f1998      0.47    0.34
year_f1999      0.66    0.34
year_f2000      0.44    0.35
year_f2001      0.37    0.32
year_f2002      0.03    0.34
year_f2003      0.57    0.34
year_f2004      0.20    0.34
year_f2005      0.66    0.39
year_f2006      1.14    0.34
year_f2007      1.40    0.33
year_f2008      0.13    0.35
year_f2009     -0.22    0.34
year_f2010     -0.06    0.33
year_f2012      0.01    0.34
year_f2013      0.17    0.33
year_f2014      0.24    0.33
year_f2015     -0.40    0.32
year_f2016     -0.01    0.35
year_f2017     -0.25    0.35
year_f2018      0.04    0.45
year_f2019      0.36    0.34
year_f2020     -0.03    0.35
year_f2021      0.02    0.38
year_f2022      0.52    0.39

Dispersion parameter: 2.82
Matérn range: 46.41
Spatial SD: 0.46
Spatiotemporal IID SD: 0.41

ML criterion at convergence: -1679.754

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameters don't look unreasonably large
Simulating ■■■■■■■■■■■                       34% | ETA:  2s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

Simulating ■■■■                               9% | ETA: 10s
Simulating ■■■■■■■                           18% | ETA:  9s
Simulating ■■■■■■■■■■■■■■                    44% | ETA:  6s
Simulating ■■■■■■■■■■■■■■■■■■■■■■            70% | ETA:  3s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■    96% | ETA:  0s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

filter: no rows removed
filter: no rows removed
mutate: new variable 'pred_dens' (double) with 135,981 unique values and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: no rows removed
filter: removed 131,292 rows (97%), 4,689 rows remaining
mutate: changed 29 values (100%) of 'type' (0 new NA)
filter: no rows removed
mutate: new variable 'species' (character) with one unique value and 0% NA
summarise: now 29 rows and 2 columns, ungrouped
Joining with `by = join_by(year)`left_join: added one column (pred_dens_sum)
           > rows only in x    0
           > rows only in y  ( 0)
           > matched rows     29
           >                 ====
           > rows total       29
mutate: changed 29 values (100%) of 'est' (0 new NA)
        changed 29 values (100%) of 'lwr' (0 new NA)
        changed 29 values (100%) of 'upr' (0 new NA)
mutate: changed 29 values (100%) of 'type' (0 new NA)
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA

[1] "Pholis gunnellus"

filter: removed 38,692 rows (88%), 5,081 rows remaining
filter: removed 16,734 rows (92%), 1,530 rows remaining

Spatial model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc + chl_sc + year_f + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: gengamma(link = 'log')
 
            coef.est coef.se
(Intercept)     2.42    0.04
temp_sc        -0.12    0.02
chl_sc          0.00    0.01
year_f1994     -0.07    0.04
year_f1995     -0.02    0.04
year_f1996     -0.23    0.06
year_f1997     -0.24    0.06
year_f1998      0.22    0.04
year_f1999      0.08    0.05
year_f2001      0.39    0.05
year_f2002     -0.02    0.04
year_f2003     -0.19    0.04
year_f2004      0.16    0.04
year_f2005      0.23    0.04
year_f2006      0.18    0.04
year_f2007      0.13    0.06
year_f2008     -0.02    0.05
year_f2009      0.09    0.04
year_f2010     -0.04    0.05
year_f2012      0.19    0.05
year_f2013      0.23    0.04
year_f2014      0.26    0.05
year_f2015      0.09    0.04
year_f2016      0.14    0.04
year_f2017      0.00    0.04
year_f2018      0.28    0.05
year_f2019      0.30    0.04
year_f2020      0.19    0.05
year_f2021      0.17    0.05
year_f2022     -0.04    0.04
yday_ct         0.01    0.00

Dispersion parameter: 0.19
Generalized gamma Q: 0.13
Matérn range: 13.50
Spatial SD: 0.14
ML criterion at convergence: 12330.722

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: no_m2 ~ year_f
Mesh: mesh_dens (isotropic covariance)
Time column: year
Data: db
Family: delta_gamma(link1 = 'log', link2 = 'log', type = 'poisson-link')

Delta/hurdle model 1: -----------------------------------
Family: binomial(link = 'log') 
            coef.est coef.se
(Intercept)    -0.20    0.92
year_f1994     -0.78    0.50
year_f1995     -0.95    0.49
year_f1996     -1.90    0.53
year_f1997     -2.05    0.57
year_f1998     -0.66    0.48
year_f1999     -1.68    0.51
year_f2000     -4.25    0.84
year_f2001     -0.66    0.46
year_f2002     -1.51    0.50
year_f2003     -1.15    0.47
year_f2004     -0.26    0.48
year_f2005     -0.16    0.51
year_f2006     -0.22    0.50
year_f2007     -1.19    0.48
year_f2008     -1.07    0.48
year_f2009     -0.06    0.50
year_f2010     -0.50    0.48
year_f2012     -1.16    0.48
year_f2013     -0.39    0.48
year_f2014     -1.44    0.48
year_f2015     -0.61    0.47
year_f2016     -1.16    0.48
year_f2017     -1.26    0.49
year_f2018     -0.61    0.51
year_f2019     -0.56    0.48
year_f2020     -0.31    0.49
year_f2021     -0.35    0.49
year_f2022     -0.98    0.52

Matérn range: 169.05
Spatial SD: 1.18
Spatiotemporal IID SD: 0.26

Delta/hurdle model 2: -----------------------------------
Family: Gamma(link = 'log') 
            coef.est coef.se
(Intercept)    -4.23    0.30
year_f1994      0.55    0.42
year_f1995      0.12    0.37
year_f1996     -0.46    0.39
year_f1997      0.41    0.42
year_f1998      0.99    0.38
year_f1999      0.51    0.38
year_f2000      0.59    0.61
year_f2001     -0.22    0.36
year_f2002     -0.02    0.38
year_f2003     -0.13    0.37
year_f2004      0.27    0.36
year_f2005      0.28    0.42
year_f2006      0.10    0.39
year_f2007      0.23    0.37
year_f2008     -0.78    0.35
year_f2009     -0.24    0.41
year_f2010     -0.32    0.36
year_f2012     -0.25    0.37
year_f2013     -0.49    0.38
year_f2014     -0.68    0.36
year_f2015     -0.91    0.37
year_f2016     -0.21    0.37
year_f2017     -0.35    0.38
year_f2018     -0.44    0.41
year_f2019     -0.30    0.38
year_f2020     -1.01    0.39
year_f2021     -0.80    0.39
year_f2022     -0.27    0.42

Dispersion parameter: 3.85
Matérn range: 39.09
Spatial SD: 0.26
Spatiotemporal IID SD: 0.47

ML criterion at convergence: -1215.117

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameters don't look unreasonably large

Simulating ■■                                 3% | ETA: 31s
Simulating ■■■■                               9% | ETA: 30s
Simulating ■■■■■■                            17% | ETA: 28s
Simulating ■■■■■■■■■                         26% | ETA: 25s
Simulating ■■■■■■■■■■■■                      35% | ETA: 22s
Simulating ■■■■■■■■■■■■■■                    44% | ETA: 19s
Simulating ■■■■■■■■■■■■■■■■■                 53% | ETA: 16s
Simulating ■■■■■■■■■■■■■■■■■■■■              62% | ETA: 13s
Simulating ■■■■■■■■■■■■■■■■■■■■■■            71% | ETA: 10s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■         80% | ETA:  7s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      88% | ETA:  4s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■    97% | ETA:  1s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

filter: removed 4,689 rows (3%), 131,292 rows remaining
filter: no rows removed
mutate: new variable 'pred_dens' (double) with 135,981 unique values and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 4,689 rows (3%), 131,292 rows remaining
filter: removed 131,292 rows (97%), 4,689 rows remaining
mutate: changed 28 values (100%) of 'type' (0 new NA)
filter: no rows removed
mutate: new variable 'species' (character) with one unique value and 0% NA
summarise: now 29 rows and 2 columns, ungrouped
Joining with `by = join_by(year)`left_join: added one column (pred_dens_sum)
           > rows only in x    0
           > rows only in y  ( 1)
           > matched rows     28
           >                 ====
           > rows total       28
mutate: changed 28 values (100%) of 'est' (0 new NA)
        changed 28 values (100%) of 'lwr' (0 new NA)
        changed 28 values (100%) of 'upr' (0 new NA)
mutate: changed 28 values (100%) of 'type' (0 new NA)
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA

[1] "Pomatoschistus sp"

filter: removed 41,851 rows (96%), 1,922 rows remaining
filter: removed 16,734 rows (92%), 1,530 rows remaining

Spatial model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc + chl_sc + year_f + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: gengamma(link = 'log')
 
            coef.est coef.se
(Intercept)     3.68    0.03
temp_sc         0.02    0.02
chl_sc         -0.01    0.01
year_f1994     -0.01    0.02
year_f1995     -0.09    0.03
year_f1996      0.01    0.05
year_f1997      0.00    0.04
year_f1998      0.06    0.03
year_f1999      0.07    0.03
year_f2000      0.07    0.03
year_f2001      0.09    0.04
year_f2002      0.01    0.03
year_f2003      0.10    0.04
year_f2004      0.10    0.03
year_f2005      0.06    0.04
year_f2006      0.07    0.04
year_f2007     -0.05    0.05
year_f2008      0.07    0.05
year_f2009      0.05    0.04
year_f2010      0.15    0.06
year_f2012     -0.01    0.04
year_f2013      0.14    0.05
year_f2014      0.03    0.05
year_f2015     -0.01    0.04
year_f2016      0.05    0.05
year_f2017     -0.01    0.04
year_f2019      0.03    0.04
year_f2022     -0.06    0.05
yday_ct         0.00    0.00

Dispersion parameter: 0.15
Generalized gamma Q: 0.42
Matérn range: 47.97
Spatial SD: 0.07
ML criterion at convergence: 6313.174

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: no_m2 ~ year_f
Mesh: mesh_dens (isotropic covariance)
Time column: year
Data: db
Family: delta_gamma(link1 = 'log', link2 = 'log', type = 'poisson-link')

Delta/hurdle model 1: -----------------------------------
Family: binomial(link = 'log') 
            coef.est coef.se
(Intercept)    -0.61    0.41
year_f1994      0.77    0.45
year_f1995      0.26    0.43
year_f1996     -1.24    0.51
year_f1997     -0.08    0.46
year_f1998      0.49    0.43
year_f1999     -0.07    0.42
year_f2000     -0.90    0.45
year_f2001     -0.49    0.42
year_f2002     -0.37    0.44
year_f2003     -0.90    0.45
year_f2004     -0.61    0.44
year_f2005     -0.42    0.46
year_f2006     -0.83    0.45
year_f2007      0.06    0.43
year_f2008     -1.14    0.47
year_f2009     -1.05    0.46
year_f2010     -1.81    0.53
year_f2012     -0.08    0.42
year_f2013     -0.86    0.45
year_f2014      0.43    0.42
year_f2015     -0.10    0.41
year_f2016     -1.25    0.47
year_f2017     -0.95    0.45
year_f2018     -1.70    0.53
year_f2019     -1.03    0.46
year_f2020     -3.30    0.82
year_f2021     -2.46    0.64
year_f2022     -1.08    0.49

Matérn range: 75.65
Spatial SD: 0.78
Spatiotemporal IID SD: 0.50

Delta/hurdle model 2: -----------------------------------
Family: Gamma(link = 'log') 
            coef.est coef.se
(Intercept)    -4.89    0.20
year_f1994      0.31    0.32
year_f1995      1.10    0.31
year_f1996      0.74    0.35
year_f1997      0.76    0.32
year_f1998      0.66    0.30
year_f1999      0.69    0.26
year_f2000      0.69    0.28
year_f2001      0.47    0.26
year_f2002      0.89    0.27
year_f2003      0.64    0.27
year_f2004      0.56    0.27
year_f2005      0.14    0.28
year_f2006      0.57    0.28
year_f2007      0.09    0.26
year_f2008      0.31    0.29
year_f2009      0.23    0.28
year_f2010      0.39    0.32
year_f2012      0.15    0.26
year_f2013      0.06    0.28
year_f2014     -0.09    0.26
year_f2015     -0.34    0.27
year_f2016      0.15    0.28
year_f2017      0.37    0.27
year_f2018     -0.02    0.33
year_f2019      0.17    0.28
year_f2020     -0.37    0.55
year_f2021      0.00    0.40
year_f2022      0.32    0.30

Dispersion parameter: 4.80
Matérn range: 14.39
Spatial SD: 0.25
Spatiotemporal IID SD: 0.51

ML criterion at convergence: -878.550

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameters don't look unreasonably large

Simulating ■■                                 3% | ETA: 31s
Simulating ■■■                                7% | ETA: 31s
Simulating ■■■■■■                            15% | ETA: 29s
Simulating ■■■■■■■■                          24% | ETA: 26s
Simulating ■■■■■■■■■■■                       33% | ETA: 23s
Simulating ■■■■■■■■■■■■■■                    42% | ETA: 20s
Simulating ■■■■■■■■■■■■■■■■                  51% | ETA: 17s
Simulating ■■■■■■■■■■■■■■■■■■■               59% | ETA: 14s
Simulating ■■■■■■■■■■■■■■■■■■■■■             68% | ETA: 11s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■          77% | ETA:  8s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■       86% | ETA:  5s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■     95% | ETA:  2s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

filter: removed 14,067 rows (10%), 121,914 rows remaining
filter: no rows removed
mutate: new variable 'pred_dens' (double) with 135,981 unique values and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 14,067 rows (10%), 121,914 rows remaining
filter: removed 131,292 rows (97%), 4,689 rows remaining
mutate: changed 26 values (100%) of 'type' (0 new NA)
filter: no rows removed
mutate: new variable 'species' (character) with one unique value and 0% NA
summarise: now 29 rows and 2 columns, ungrouped
Joining with `by = join_by(year)`left_join: added one column (pred_dens_sum)
           > rows only in x    0
           > rows only in y  ( 3)
           > matched rows     26
           >                 ====
           > rows total       26
mutate: changed 26 values (100%) of 'est' (0 new NA)
        changed 26 values (100%) of 'lwr' (0 new NA)
        changed 26 values (100%) of 'upr' (0 new NA)
mutate: changed 26 values (100%) of 'type' (0 new NA)
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA

[1] "Ammodytidae"

filter: removed 42,616 rows (97%), 1,157 rows remaining
filter: removed 16,734 rows (92%), 1,530 rows remaining

Spatial model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc + chl_sc + year_f + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: gengamma(link = 'log')
 
            coef.est coef.se
(Intercept)     3.62    0.29
temp_sc         0.00    0.04
chl_sc         -0.03    0.02
year_f1994     -0.06    0.07
year_f1995      0.13    0.10
year_f1997     -0.21    0.12
year_f1998     -0.10    0.08
year_f1999      0.01    0.08
year_f2000     -0.36    0.08
year_f2001      0.06    0.10
year_f2002     -0.16    0.08
year_f2003     -0.23    0.10
year_f2004     -0.14    0.07
year_f2005     -0.07    0.09
year_f2006     -0.08    0.08
year_f2007     -0.04    0.12
year_f2008     -0.09    0.08
year_f2009     -0.01    0.09
year_f2010     -0.24    0.12
year_f2012     -0.21    0.10
year_f2013     -0.09    0.09
year_f2014     -0.12    0.10
year_f2015     -0.12    0.10
year_f2016     -0.07    0.10
year_f2017     -0.24    0.07
year_f2018     -0.09    0.09
year_f2019     -0.14    0.08
year_f2020     -0.14    0.10
year_f2021     -0.25    0.09
year_f2022     -0.23    0.09
yday_ct         0.00    0.00

Dispersion parameter: 0.15
Generalized gamma Q: 2.68
Matérn range: 161.91
Spatial SD: 0.41
ML criterion at convergence: 4220.239

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: no_m2 ~ year_f
Mesh: mesh_dens (isotropic covariance)
Time column: year
Data: db
Family: delta_gamma(link1 = 'log', link2 = 'log', type = 'poisson-link')

Delta/hurdle model 1: -----------------------------------
Family: binomial(link = 'log') 
            coef.est coef.se
(Intercept)    -0.63    0.41
year_f1994      0.08    0.49
year_f1995     -0.66    0.51
year_f1996     -1.94    0.67
year_f1997     -0.66    0.53
year_f1998     -0.73    0.51
year_f1999     -0.83    0.51
year_f2000     -0.55    0.50
year_f2001     -1.33    0.53
year_f2002     -0.28    0.49
year_f2003     -0.78    0.50
year_f2004     -0.19    0.47
year_f2005     -0.96    0.56
year_f2006     -1.05    0.53
year_f2007     -0.92    0.52
year_f2008      0.58    0.47
year_f2009     -0.61    0.50
year_f2010     -2.05    0.66
year_f2012     -0.49    0.49
year_f2013     -0.03    0.47
year_f2014      0.11    0.46
year_f2015     -0.63    0.48
year_f2016     -0.56    0.48
year_f2017      0.62    0.46
year_f2018      0.34    0.47
year_f2019     -0.19    0.47
year_f2020      1.16    0.47
year_f2021      0.18    0.49
year_f2022      0.19    0.48

Matérn range: 112.35
Spatial SD: 0.47
Spatiotemporal IID SD: 0.45

Delta/hurdle model 2: -----------------------------------
Family: Gamma(link = 'log') 
            coef.est coef.se
(Intercept)    -4.81    0.20
year_f1994      0.10    0.24
year_f1995      0.57    0.26
year_f1996      0.15    0.35
year_f1997      0.69    0.27
year_f1998      0.38    0.25
year_f1999      0.58    0.25
year_f2000      0.64    0.26
year_f2001      0.35    0.26
year_f2002      0.73    0.24
year_f2003      0.25    0.25
year_f2004      0.21    0.23
year_f2005      0.71    0.30
year_f2006      0.43    0.27
year_f2007      0.17    0.26
year_f2008      0.17    0.23
year_f2009     -0.06    0.24
year_f2010      0.04    0.35
year_f2012      0.06    0.25
year_f2013     -0.03    0.23
year_f2014      0.11    0.23
year_f2015     -0.27    0.24
year_f2016     -0.01    0.23
year_f2017      0.65    0.26
year_f2018      0.02    0.25
year_f2019      0.00    0.23
year_f2020     -0.62    0.25
year_f2021      0.18    0.25
year_f2022     -0.08    0.24

Dispersion parameter: 3.69
Matérn range: 62.56
Spatial SD: 0.33
Spatiotemporal IID SD: 0.22

ML criterion at convergence: -920.031

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameters don't look unreasonably large

Simulating ■■■■                               9% | ETA: 10s
Simulating ■■■■■                             15% | ETA: 10s
Simulating ■■■■■■■■■■■■■                     42% | ETA:  7s
Simulating ■■■■■■■■■■■■■■■■■■■■■■            69% | ETA:  4s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■    95% | ETA:  1s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

filter: removed 4,689 rows (3%), 131,292 rows remaining
filter: no rows removed
mutate: new variable 'pred_dens' (double) with 135,981 unique values and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 4,689 rows (3%), 131,292 rows remaining
filter: removed 131,292 rows (97%), 4,689 rows remaining
mutate: changed 28 values (100%) of 'type' (0 new NA)
filter: no rows removed
mutate: new variable 'species' (character) with one unique value and 0% NA
summarise: now 29 rows and 2 columns, ungrouped
Joining with `by = join_by(year)`left_join: added one column (pred_dens_sum)
           > rows only in x    0
           > rows only in y  ( 1)
           > matched rows     28
           >                 ====
           > rows total       28
mutate: changed 28 values (100%) of 'est' (0 new NA)
        changed 28 values (100%) of 'lwr' (0 new NA)
        changed 28 values (100%) of 'upr' (0 new NA)
mutate: changed 28 values (100%) of 'type' (0 new NA)
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA

[1] "Chirolophis ascanii"

filter: removed 42,719 rows (98%), 1,054 rows remaining
filter: removed 16,734 rows (92%), 1,530 rows remaining

Spatial model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc + chl_sc + year_f + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: gengamma(link = 'log')
 
            coef.est coef.se
(Intercept)     2.74    0.07
temp_sc        -0.07    0.03
chl_sc          0.04    0.02
year_f1994     -0.10    0.07
year_f1995      0.13    0.07
year_f1996     -0.40    0.09
year_f1997     -0.27    0.11
year_f1998      0.10    0.07
year_f1999     -0.03    0.07
year_f2001      0.22    0.09
year_f2002     -0.22    0.08
year_f2003     -0.22    0.09
year_f2004     -0.01    0.07
year_f2005      0.10    0.09
year_f2006     -0.05    0.08
year_f2007      0.06    0.11
year_f2008     -0.05    0.10
year_f2009     -0.10    0.07
year_f2010     -0.19    0.09
year_f2012      0.01    0.10
year_f2013      0.03    0.08
year_f2014     -0.04    0.11
year_f2015      0.01    0.09
year_f2016     -0.09    0.10
year_f2017     -0.20    0.08
year_f2018     -0.03    0.09
year_f2019     -0.17    0.08
year_f2020     -0.02    0.10
year_f2021     -0.08    0.08
year_f2022     -0.19    0.08
yday_ct         0.00    0.00

Dispersion parameter: 0.24
Generalized gamma Q: 0.41
Matérn range: 9.62
Spatial SD: 0.12
ML criterion at convergence: 2970.371

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: no_m2 ~ year_f
Mesh: mesh_dens (isotropic covariance)
Time column: year
Data: db
Family: delta_gamma(link1 = 'log', link2 = 'log', type = 'poisson-link')

Delta/hurdle model 1: -----------------------------------
Family: binomial(link = 'log') 
            coef.est coef.se
(Intercept)    -1.48    1.26
year_f1994     -0.28    0.37
year_f1995     -0.55    0.37
year_f1996     -0.58    0.36
year_f1997     -0.82    0.43
year_f1998      0.28    0.31
year_f1999     -0.25    0.33
year_f2000     -1.45    0.45
year_f2001     -0.29    0.31
year_f2002     -0.45    0.36
year_f2003     -0.85    0.38
year_f2004     -0.23    0.33
year_f2005     -1.07    0.47
year_f2006     -0.38    0.36
year_f2007     -0.91    0.40
year_f2008     -0.69    0.38
year_f2009      0.68    0.32
year_f2010     -0.20    0.34
year_f2012      0.25    0.32
year_f2013      0.01    0.32
year_f2014     -0.68    0.37
year_f2015     -0.33    0.34
year_f2016     -0.99    0.40
year_f2017     -0.31    0.35
year_f2018      0.16    0.35
year_f2019     -0.10    0.34
year_f2020      0.06    0.35
year_f2021      0.04    0.35
year_f2022      0.56    0.35

Matérn range: 254.26
Spatial SD: 1.15
Spatiotemporal IID SD: 0.00

Delta/hurdle model 2: -----------------------------------
Family: Gamma(link = 'log') 
            coef.est coef.se
(Intercept)    -4.74    0.19
year_f1994      0.05    0.29
year_f1995      0.51    0.28
year_f1996      0.03    0.28
year_f1997      0.30    0.32
year_f1998      0.56    0.26
year_f1999      0.60    0.26
year_f2000     -0.08    0.31
year_f2001      0.23    0.24
year_f2002     -0.02    0.28
year_f2003      0.33    0.28
year_f2004      0.18    0.25
year_f2005      0.44    0.33
year_f2006      0.10    0.27
year_f2007      0.07    0.29
year_f2008     -0.65    0.28
year_f2009     -0.18    0.27
year_f2010      0.06    0.26
year_f2012     -0.33    0.27
year_f2013      0.08    0.25
year_f2014     -0.50    0.27
year_f2015     -0.39    0.27
year_f2016     -0.23    0.28
year_f2017      0.19    0.27
year_f2018      0.10    0.28
year_f2019      0.12    0.26
year_f2020     -0.53    0.27
year_f2021      0.38    0.27
year_f2022     -0.18    0.30

Dispersion parameter: 6.08
Matérn range: 26.27
Spatial SD: 0.30
Spatiotemporal IID SD: 0.44

ML criterion at convergence: -877.930

See ?tidy.sdmTMB to extract these values as a data frame.

**Possible issues detected! Check output of sanity().**

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✖ `ln_tau_E` standard error may be large
ℹ `ln_tau_E` is an internal parameter affecting `sigma_E`
ℹ `sigma_E` is the spatiotemporal standard deviation
ℹ Try simplifying the model, adjusting the mesh, or adding priors

✖ `sigma_E` is smaller than 0.01
ℹ Consider omitting this part of the model

✔ Range parameters don't look unreasonably large

Simulating ■■                                 3% | ETA: 31s
Simulating ■■■■                               9% | ETA: 30s
Simulating ■■■■■■                            18% | ETA: 27s
Simulating ■■■■■■■■■                         27% | ETA: 25s
Simulating ■■■■■■■■■■■■                      36% | ETA: 22s
Simulating ■■■■■■■■■■■■■■                    45% | ETA: 19s
Simulating ■■■■■■■■■■■■■■■■■                 53% | ETA: 16s
Simulating ■■■■■■■■■■■■■■■■■■■■              62% | ETA: 13s
Simulating ■■■■■■■■■■■■■■■■■■■■■■            71% | ETA: 10s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■         80% | ETA:  7s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      88% | ETA:  4s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■    97% | ETA:  1s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

filter: removed 4,689 rows (3%), 131,292 rows remaining
filter: no rows removed
mutate: new variable 'pred_dens' (double) with 135,981 unique values and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 4,689 rows (3%), 131,292 rows remaining
filter: removed 131,292 rows (97%), 4,689 rows remaining
mutate: changed 28 values (100%) of 'type' (0 new NA)
filter: no rows removed
mutate: new variable 'species' (character) with one unique value and 0% NA
summarise: now 29 rows and 2 columns, ungrouped
Joining with `by = join_by(year)`left_join: added one column (pred_dens_sum)
           > rows only in x    0
           > rows only in y  ( 1)
           > matched rows     28
           >                 ====
           > rows total       28
mutate: changed 28 values (100%) of 'est' (0 new NA)
        changed 28 values (100%) of 'lwr' (0 new NA)
        changed 28 values (100%) of 'upr' (0 new NA)
mutate: changed 28 values (100%) of 'type' (0 new NA)
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA

[1] "Agonus cataphractus"

filter: removed 43,435 rows (99%), 338 rows remaining
filter: removed 16,850 rows (92%), 1,414 rows remaining

Warning in sqrt(diag(cov)): NaNs produced

Warning: The model may not have converged: non-positive-definite Hessian
matrix.

Spatial model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc + chl_sc + year_f + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: gengamma(link = 'log')
 
            coef.est coef.se
(Intercept)     1.90    0.09
temp_sc        -0.01    0.05
chl_sc          0.03    0.02
year_f1994      0.05    0.05
year_f1998      0.40    0.09
year_f2003      0.09    0.09
year_f2004      0.43    0.09
year_f2005      0.20    0.09
year_f2006      0.24    0.07
year_f2007      0.32    0.13
year_f2009      0.22    0.07
year_f2013      0.26    0.10
year_f2015      0.26    0.10
year_f2017      0.24    0.09
year_f2019      0.25    0.07
year_f2022      0.21    0.10
yday_ct         0.01    0.00

Dispersion parameter: 0.14
Generalized gamma Q: -0.41
Matérn range: 6.00
Spatial SD: 0.08
ML criterion at convergence: 573.864

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: no_m2 ~ year_f
Mesh: mesh_dens (isotropic covariance)
Time column: year
Data: db
Family: delta_gamma(link1 = 'log', link2 = 'log', type = 'poisson-link')

Delta/hurdle model 1: -----------------------------------

Warning in sqrt(diag(object$cov.fixed)): NaNs produced

Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced

Family: binomial(link = 'log') 
            coef.est coef.se
(Intercept)    -3.14    2.39
year_f1994      1.04    0.43
year_f1995     -2.61    1.05
year_f1996     -0.87    0.55
year_f1997     -1.52    0.78
year_f1998     -0.25    0.46
year_f1999     -1.87    0.78
year_f2001     -2.24    0.78
year_f2002     -2.17    1.05
year_f2003      0.16    0.46
year_f2004     -0.65    0.49
year_f2005      0.44    0.50
year_f2006     -0.70    0.53
year_f2007     -0.52    0.48
year_f2008     -1.37    0.62
year_f2009      0.14    0.45
year_f2010     -0.84    0.53
year_f2013      0.05    0.45
year_f2014     -1.27    0.56
year_f2015     -0.24    0.46
year_f2016     -0.97    0.56
year_f2017     -0.65    0.48
year_f2018     -2.14    1.05
year_f2019     -0.08    0.49
year_f2020     -0.47    0.53
year_f2021     -1.48    0.78
year_f2022     -0.71    0.59

Matérn range: 249.95
Spatial SD: 2.21
Spatiotemporal IID SD: 0.00

Delta/hurdle model 2: -----------------------------------

Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced

Family: Gamma(link = 'log') 
            coef.est coef.se
(Intercept)    -4.78    0.20
year_f1994     -0.24    0.30
year_f1995      1.21    0.58
year_f1996      0.26    0.31
year_f1997      0.62    0.42
year_f1998      0.32    0.28
year_f1999      1.39    0.44
year_f2001     -0.29    0.43
year_f2002      0.42    0.56
year_f2003      0.41    0.28
year_f2004      0.08    0.29
year_f2005      0.04    0.36
year_f2006      0.33    0.30
year_f2007     -0.43    0.28
year_f2008     -0.63    0.33
year_f2009     -0.01    0.28
year_f2010     -0.16    0.29
year_f2013     -0.03    0.28
year_f2014     -0.64    0.31
year_f2015     -0.61    0.29
year_f2016     -0.67    0.33
year_f2017     -0.37    0.27
year_f2018      0.32    0.56
year_f2019      0.16    0.31
year_f2020     -0.54    0.30
year_f2021     -0.47    0.43
year_f2022     -0.12    0.41

Dispersion parameter: 3.57
Matérn range: 0.00
Spatial SD: 0.00
Spatiotemporal IID SD: 0.00

ML criterion at convergence: -298.609

See ?tidy.sdmTMB to extract these values as a data frame.

**Possible issues detected! Check output of sanity().**

✔ Non-linear minimizer suggests successful convergence
✖ Non-positive-definite Hessian matrix: model may not have converged
ℹ Try simplifying the model, adjusting the mesh, or adding priors

✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001

Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced

✔ No fixed-effect standard errors are NA

Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced

✖ `ln_tau_E` standard error may be large
ℹ `ln_tau_E` is an internal parameter affecting `sigma_E`
ℹ `sigma_E` is the spatiotemporal standard deviation
ℹ Try simplifying the model, adjusting the mesh, or adding priors

Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced

✖ `sigma_E` is smaller than 0.01
ℹ Consider omitting this part of the model

✖ `sigma_O` is smaller than 0.01
ℹ Consider omitting this part of the model

✖ `sigma_E` is smaller than 0.01
ℹ Consider omitting this part of the model

✔ Range parameters don't look unreasonably large

Simulating ■■                                 3% | ETA: 29s
Simulating ■■■■■                             12% | ETA: 27s
Simulating ■■■■■■■                           22% | ETA: 25s
Simulating ■■■■■■■■■■                        31% | ETA: 21s
Simulating ■■■■■■■■■■■■■                     41% | ETA: 18s
Simulating ■■■■■■■■■■■■■■■■                  51% | ETA: 15s
Simulating ■■■■■■■■■■■■■■■■■■■               60% | ETA: 12s
Simulating ■■■■■■■■■■■■■■■■■■■■■■            70% | ETA:  9s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■         79% | ETA:  7s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      89% | ETA:  4s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■   98% | ETA:  1s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

filter: removed 70,335 rows (52%), 65,646 rows remaining
filter: removed 9,378 rows (7%), 126,603 rows remaining
mutate: new variable 'pred_dens' (double) with 126,603 unique values and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 60,957 rows (48%), 65,646 rows remaining
filter: removed 131,292 rows (97%), 4,689 rows remaining
mutate: changed 14 values (100%) of 'type' (0 new NA)
filter: removed 9,378 rows (7%), 126,603 rows remaining

Warning in sqrt(diag(cov)): NaNs produced

Warning: The model may not have converged: non-positive-definite Hessian
matrix.

mutate: new variable 'species' (character) with one unique value and 0% NA
summarise: now 27 rows and 2 columns, ungrouped
Joining with `by = join_by(year)`left_join: added one column (pred_dens_sum)
           > rows only in x    0
           > rows only in y  (13)
           > matched rows     14
           >                 ====
           > rows total       14
mutate: changed 14 values (100%) of 'est' (0 new NA)
        changed 14 values (100%) of 'lwr' (0 new NA)
        changed 14 values (100%) of 'upr' (0 new NA)
mutate: changed 14 values (100%) of 'type' (0 new NA)
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA

Warning in sqrt(diag(object$cov.fixed)): NaNs produced

Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced

filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA

[1] "Microstomus kitt"

filter: removed 43,164 rows (99%), 609 rows remaining
filter: removed 16,790 rows (92%), 1,474 rows remaining

Spatial model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc + chl_sc + year_f + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: gengamma(link = 'log')
 
            coef.est coef.se
(Intercept)     3.21    0.10
temp_sc        -0.06    0.04
chl_sc          0.05    0.03
year_f1994      0.08    0.10
year_f1995     -0.08    0.08
year_f1998      0.11    0.11
year_f1999     -0.05    0.10
year_f2000      0.18    0.12
year_f2001      0.11    0.13
year_f2002      0.17    0.11
year_f2003      0.01    0.13
year_f2006      0.20    0.11
year_f2008      0.15    0.12
year_f2009      0.00    0.12
year_f2012     -0.09    0.14
year_f2014      0.02    0.15
year_f2015     -0.07    0.12
year_f2016      0.13    0.14
year_f2017     -0.03    0.11
year_f2018     -0.22    0.12
year_f2019     -0.21    0.11
year_f2020     -0.16    0.13
year_f2021     -0.58    0.11
year_f2022     -0.06    0.13
yday_ct         0.01    0.00

Dispersion parameter: 0.25
Generalized gamma Q: 0.18
Matérn range: 24.61
Spatial SD: 0.07
ML criterion at convergence: 1880.858

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: no_m2 ~ year_f
Mesh: mesh_dens (isotropic covariance)
Time column: year
Data: db
Family: delta_gamma(link1 = 'log', link2 = 'log', type = 'poisson-link')

Delta/hurdle model 1: -----------------------------------
Family: binomial(link = 'log') 
            coef.est coef.se
(Intercept)    -2.03    0.45
year_f1994      0.90    0.52
year_f1995      1.47    0.49
year_f1996     -2.12    1.10
year_f1997     -0.53    0.64
year_f1998      0.01    0.54
year_f1999      0.55    0.50
year_f2000     -0.19    0.55
year_f2001     -0.46    0.57
year_f2002     -0.01    0.54
year_f2003     -0.49    0.57
year_f2004     -1.05    0.63
year_f2005     -0.66    0.65
year_f2006      0.02    0.54
year_f2008      0.37    0.53
year_f2009     -0.51    0.58
year_f2010     -2.52    1.10
year_f2012      0.74    0.50
year_f2013     -0.91    0.63
year_f2014      0.81    0.49
year_f2015      0.12    0.54
year_f2016     -0.17    0.54
year_f2017      1.19    0.48
year_f2018     -1.11    0.68
year_f2019      0.72    0.51
year_f2020     -0.36    0.62
year_f2021      1.02    0.51
year_f2022      0.32    0.56

Matérn range: 47.99
Spatial SD: 0.85
Spatiotemporal IID SD: 0.64

Delta/hurdle model 2: -----------------------------------
Family: Gamma(link = 'log') 
            coef.est coef.se
(Intercept)    -4.88    0.20
year_f1994      0.15    0.27
year_f1995      0.23    0.26
year_f1996     -0.02    0.52
year_f1997      0.26    0.29
year_f1998      0.16    0.24
year_f1999      0.73    0.23
year_f2000      0.35    0.25
year_f2001      0.68    0.25
year_f2002      0.37    0.24
year_f2003      0.55    0.24
year_f2004     -0.11    0.28
year_f2005      0.31    0.29
year_f2006      0.25    0.26
year_f2008      0.26    0.26
year_f2009      0.40    0.26
year_f2010      0.05    0.54
year_f2012     -0.08    0.22
year_f2013     -0.02    0.28
year_f2014     -0.25    0.24
year_f2015     -0.20    0.24
year_f2016     -0.07    0.24
year_f2017      0.01    0.23
year_f2018      1.53    0.32
year_f2019      0.13    0.26
year_f2020     -0.18    0.28
year_f2021     -0.05    0.26
year_f2022      0.17    0.25

Dispersion parameter: 6.95
Matérn range: 25.24
Spatial SD: 0.37
Spatiotemporal IID SD: 0.00

ML criterion at convergence: -441.924

See ?tidy.sdmTMB to extract these values as a data frame.

**Possible issues detected! Check output of sanity().**

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✖ `ln_tau_E` standard error may be large
ℹ `ln_tau_E` is an internal parameter affecting `sigma_E`
ℹ `sigma_E` is the spatiotemporal standard deviation
ℹ Try simplifying the model, adjusting the mesh, or adding priors

✖ `sigma_E` is smaller than 0.01
ℹ Consider omitting this part of the model

✔ Range parameters don't look unreasonably large

Simulating ■■                                 3% | ETA: 31s
Simulating ■■■■                               9% | ETA: 30s
Simulating ■■■■■■                            18% | ETA: 27s
Simulating ■■■■■■■■■                         27% | ETA: 25s
Simulating ■■■■■■■■■■■■                      35% | ETA: 22s
Simulating ■■■■■■■■■■■■■■                    44% | ETA: 19s
Simulating ■■■■■■■■■■■■■■■■■                 53% | ETA: 16s
Simulating ■■■■■■■■■■■■■■■■■■■■              62% | ETA: 13s
Simulating ■■■■■■■■■■■■■■■■■■■■■■            71% | ETA: 10s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■         80% | ETA:  7s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      89% | ETA:  4s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■    98% | ETA:  1s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

filter: removed 32,823 rows (24%), 103,158 rows remaining
filter: removed 4,689 rows (3%), 131,292 rows remaining
mutate: new variable 'pred_dens' (double) with 131,292 unique values and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 28,134 rows (21%), 103,158 rows remaining
filter: removed 131,292 rows (97%), 4,689 rows remaining
mutate: changed 22 values (100%) of 'type' (0 new NA)
filter: removed 4,689 rows (3%), 131,292 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
summarise: now 28 rows and 2 columns, ungrouped
Joining with `by = join_by(year)`left_join: added one column (pred_dens_sum)
           > rows only in x    0
           > rows only in y  ( 6)
           > matched rows     22
           >                 ====
           > rows total       22
mutate: changed 22 values (100%) of 'est' (0 new NA)
        changed 22 values (100%) of 'lwr' (0 new NA)
        changed 22 values (100%) of 'upr' (0 new NA)
mutate: changed 22 values (100%) of 'type' (0 new NA)
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA

[1] "Syngnathus rostellatus"

filter: removed 42,730 rows (98%), 1,043 rows remaining
filter: removed 16,734 rows (92%), 1,530 rows remaining

Spatial model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc + chl_sc + year_f + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: gengamma(link = 'log')
 
            coef.est coef.se
(Intercept)     4.24    0.08
temp_sc        -0.02    0.02
chl_sc         -0.06    0.02
year_f1998     -0.02    0.09
year_f1999     -0.03    0.08
year_f2000     -0.10    0.09
year_f2001      0.06    0.10
year_f2003     -0.18    0.08
year_f2004     -0.17    0.08
year_f2005     -0.01    0.09
year_f2006     -0.12    0.09
year_f2007      0.10    0.10
year_f2008     -0.14    0.08
year_f2009     -0.13    0.09
year_f2012      0.02    0.09
year_f2013      0.06    0.09
year_f2014      0.08    0.09
year_f2015      0.07    0.08
year_f2016     -0.07    0.08
year_f2018     -0.11    0.09
year_f2019     -0.14    0.09
year_f2020     -0.26    0.09
year_f2021     -0.05    0.09
year_f2022     -0.33    0.09
yday_ct         0.01    0.00

Dispersion parameter: 0.22
Generalized gamma Q: 0.15
Matérn range: 22.39
Spatial SD: 0.04
ML criterion at convergence: 4247.354

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: no_m2 ~ year_f
Mesh: mesh_dens (isotropic covariance)
Time column: year
Data: db
Family: delta_gamma(link1 = 'log', link2 = 'log', type = 'poisson-link')

Delta/hurdle model 1: -----------------------------------
Family: binomial(link = 'log') 
            coef.est coef.se
(Intercept)    -2.38    0.54
year_f1994      0.13    0.60
year_f1995     -0.91    0.74
year_f1996     -1.73    1.10
year_f1997     -1.88    1.10
year_f1998      0.62    0.53
year_f1999      0.96    0.52
year_f2000      0.55    0.53
year_f2001      0.28    0.53
year_f2002     -0.43    0.60
year_f2003      1.70    0.50
year_f2004      2.45    0.50
year_f2005      1.54    0.54
year_f2006      0.83    0.53
year_f2007      1.11    0.51
year_f2008      1.53    0.51
year_f2009      0.30    0.56
year_f2010     -0.25    0.61
year_f2012      1.44    0.51
year_f2013      1.32    0.52
year_f2014      1.45    0.50
year_f2015      2.07    0.49
year_f2016      1.46    0.50
year_f2017     -1.59    0.84
year_f2018      1.01    0.53
year_f2019      0.54    0.54
year_f2020      1.56    0.52
year_f2021      0.81    0.54
year_f2022      1.24    0.54

Matérn range: 73.63
Spatial SD: 0.95
Spatiotemporal IID SD: 0.40

Delta/hurdle model 2: -----------------------------------
Family: Gamma(link = 'log') 
            coef.est coef.se
(Intercept)    -5.01    0.29
year_f1994      0.33    0.40
year_f1995      0.74    0.51
year_f1996      0.67    0.75
year_f1997      0.62    0.68
year_f1998      0.77    0.34
year_f1999      0.85    0.34
year_f2000      0.53    0.34
year_f2001      0.89    0.34
year_f2002      0.34    0.39
year_f2003      1.09    0.34
year_f2004      0.32    0.35
year_f2005      0.35    0.35
year_f2006      0.25    0.34
year_f2007      0.24    0.33
year_f2008      0.34    0.34
year_f2009      0.45    0.36
year_f2010      0.31    0.38
year_f2012      0.01    0.34
year_f2013      0.56    0.34
year_f2014      0.37    0.34
year_f2015     -0.04    0.33
year_f2016      0.44    0.34
year_f2017      0.98    0.52
year_f2018      0.45    0.35
year_f2019      0.26    0.35
year_f2020     -0.27    0.34
year_f2021      0.13    0.35
year_f2022      0.09    0.35

Dispersion parameter: 5.41
Matérn range: 20.04
Spatial SD: 0.31
Spatiotemporal IID SD: 0.52

ML criterion at convergence: -716.474

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameters don't look unreasonably large

Simulating ■■■■                               9% | ETA: 10s
Simulating ■■■■■■■■■                         28% | ETA:  8s
Simulating ■■■■■■■■■■■■■■■■■■                55% | ETA:  5s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■        82% | ETA:  2s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

filter: removed 32,823 rows (24%), 103,158 rows remaining
filter: no rows removed
mutate: new variable 'pred_dens' (double) with 135,981 unique values and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 32,823 rows (24%), 103,158 rows remaining
filter: removed 131,292 rows (97%), 4,689 rows remaining
mutate: changed 22 values (100%) of 'type' (0 new NA)
filter: no rows removed
mutate: new variable 'species' (character) with one unique value and 0% NA
summarise: now 29 rows and 2 columns, ungrouped
Joining with `by = join_by(year)`left_join: added one column (pred_dens_sum)
           > rows only in x    0
           > rows only in y  ( 7)
           > matched rows     22
           >                 ====
           > rows total       22
mutate: changed 22 values (100%) of 'est' (0 new NA)
        changed 22 values (100%) of 'lwr' (0 new NA)
        changed 22 values (100%) of 'upr' (0 new NA)
mutate: changed 22 values (100%) of 'type' (0 new NA)
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA

[1] "Myoxocephalus scorpius"

filter: removed 43,244 rows (99%), 529 rows remaining
filter: removed 16,936 rows (93%), 1,328 rows remaining

Warning in sqrt(diag(cov)): NaNs produced

Warning: The model may not have converged: non-positive-definite Hessian
matrix.

Spatial model fit by ML ['sdmTMB']
Formula: length_mm ~ temp_sc + chl_sc + year_f + yday_ct
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: gengamma(link = 'log')
 
            coef.est coef.se
(Intercept)     1.88    0.07
temp_sc        -0.13    0.03
chl_sc         -0.02    0.02
year_f1999      0.17    0.06
year_f2003     -0.05    0.07
year_f2004      0.22    0.07
year_f2005      0.27    0.06
year_f2006      0.11    0.06
year_f2007      0.40    0.09
year_f2008      0.12    0.08
year_f2009      0.16    0.05
year_f2010     -0.02    0.07
year_f2013      0.23    0.07
year_f2015      0.10    0.09
year_f2016      0.18    0.08
year_f2017      0.13    0.07
year_f2019      0.25    0.06
year_f2020      0.14    0.08
year_f2022     -0.03    0.07
yday_ct         0.01    0.00

Dispersion parameter: 0.14
Generalized gamma Q: -0.12
Matérn range: 32.25
Spatial SD: 0.05
ML criterion at convergence: 859.102

See ?tidy.sdmTMB to extract these values as a data frame.

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: no_m2 ~ year_f
Mesh: mesh_dens (isotropic covariance)
Time column: year
Data: db
Family: delta_gamma(link1 = 'log', link2 = 'log', type = 'poisson-link')

Delta/hurdle model 1: -----------------------------------

Warning in sqrt(diag(object$cov.fixed)): NaNs produced

Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced

Family: binomial(link = 'log') 
            coef.est coef.se
(Intercept)    -2.25    1.92
year_f1995     -1.75    0.81
year_f1996     -1.32    0.68
year_f1998     -0.74    0.55
year_f1999     -0.05    0.49
year_f2001     -2.77    1.07
year_f2003      0.26    0.47
year_f2004     -0.45    0.51
year_f2005      0.95    0.46
year_f2006     -0.07    0.50
year_f2007     -0.25    0.50
year_f2008     -0.61    0.55
year_f2009      0.95    0.44
year_f2010     -0.09    0.50
year_f2012     -0.99    0.58
year_f2013      0.35    0.46
year_f2014     -1.48    0.69
year_f2015     -1.55    0.68
year_f2016     -0.55    0.53
year_f2017      0.33    0.45
year_f2018     -0.53    0.57
year_f2019      1.35    0.43
year_f2020      0.95    0.46
year_f2021     -1.47    0.80
year_f2022     -0.41    0.58

Matérn range: 362.39
Spatial SD: 1.23
Spatiotemporal IID SD: 0.00

Delta/hurdle model 2: -----------------------------------

Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced

Family: Gamma(link = 'log') 
            coef.est coef.se
(Intercept)    -4.86    0.28
year_f1995      0.58    0.50
year_f1996     -0.61    0.41
year_f1998      0.32    0.37
year_f1999      0.49    0.35
year_f2001     -0.82    0.66
year_f2003      0.03    0.34
year_f2004      0.16    0.36
year_f2005      0.49    0.34
year_f2006      0.18    0.34
year_f2007      0.64    0.35
year_f2008     -0.39    0.36
year_f2009     -0.22    0.33
year_f2010     -0.15    0.34
year_f2012     -0.51    0.38
year_f2013     -0.22    0.35
year_f2014     -0.31    0.43
year_f2015     -0.25    0.44
year_f2016     -0.25    0.36
year_f2017     -0.29    0.33
year_f2018      0.22    0.37
year_f2019     -0.22    0.34
year_f2020     -0.41    0.34
year_f2021      0.48    0.52
year_f2022      0.41    0.42

Dispersion parameter: 2339999.79
Matérn range: 21.50
Spatial SD: 0.31
Spatiotemporal IID SD: 0.53

ML criterion at convergence: -361.113

See ?tidy.sdmTMB to extract these values as a data frame.

**Possible issues detected! Check output of sanity().**

✖ Non-linear minimizer did not converge: do not trust this model!
ℹ Try simplifying the model, adjusting the mesh, or adding priors

✖ Non-positive-definite Hessian matrix: model may not have converged
ℹ Try simplifying the model, adjusting the mesh, or adding priors

✔ No extreme or very small eigenvalues detected
✖ `b_j` gradient > 0.001
ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the model

✖ `b_j` gradient > 0.001
ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the model

✖ `b_j2` gradient > 0.001
ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the model

✖ `b_j2` gradient > 0.001
ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the model

✖ `ln_tau_O` gradient > 0.001
ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the model

✖ `ln_tau_E` gradient > 0.001
ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the model

✖ `ln_kappa` gradient > 0.001
ℹ See ?run_extra_optimization(), standardize covariates, and/or simplify the model

Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced

✔ No fixed-effect standard errors are NA

Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced

✔ No standard errors look unreasonably large

Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced

✖ `sigma_E` is smaller than 0.01
ℹ Consider omitting this part of the model

✔ Range parameters don't look unreasonably large

Simulating ■■                                 4% | ETA: 26s
Simulating ■■■                                6% | ETA: 26s
Simulating ■■■■■■                            16% | ETA: 24s
Simulating ■■■■■■■■■                         26% | ETA: 21s
Simulating ■■■■■■■■■■■■                      37% | ETA: 18s
Simulating ■■■■■■■■■■■■■■■                   48% | ETA: 15s
Simulating ■■■■■■■■■■■■■■■■■■                58% | ETA: 12s
Simulating ■■■■■■■■■■■■■■■■■■■■■■            69% | ETA:  9s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■         80% | ETA:  6s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      90% | ETA:  3s
Simulating ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s

filter: removed 56,268 rows (41%), 79,713 rows remaining
filter: removed 18,756 rows (14%), 117,225 rows remaining
mutate: new variable 'pred_dens' (double) with 117,225 unique values and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 37,512 rows (32%), 79,713 rows remaining
filter: removed 131,292 rows (97%), 4,689 rows remaining
mutate: changed 17 values (100%) of 'type' (0 new NA)
filter: removed 18,756 rows (14%), 117,225 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
summarise: now 25 rows and 2 columns, ungrouped
Joining with `by = join_by(year)`left_join: added one column (pred_dens_sum)
           > rows only in x    0
           > rows only in y  ( 8)
           > matched rows     17
           >                 ====
           > rows total       17
mutate: changed 17 values (100%) of 'est' (0 new NA)
        changed 17 values (100%) of 'lwr' (0 new NA)
        changed 17 values (100%) of 'upr' (0 new NA)
mutate: changed 17 values (100%) of 'type' (0 new NA)
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA

Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced
Warning in sqrt(diag(object$cov.fixed)): NaNs produced

filter: removed 2 rows (67%), one row remaining
mutate: new variable 'cutoff' (double) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 2 rows (67%), one row remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA

pars <- bind_rows(pars_list)
res <- bind_rows(res_list)
res_dens <- bind_rows(res_dens_list)
spatial <- bind_rows(spatial_list)
spatial_density <- bind_rows(spatial_density_list)
range <- bind_rows(range_list)
genq <- bind_rows(genq_list)
index <- bind_rows(index_list)
dens_index <- bind_rows(dens_index_list)

write_csv(pars, paste0(home, "/output/pars.csv"))
write_csv(res, paste0(home, "/output/res.csv"))
write_csv(res_dens, paste0(home, "/output/res_dens.csv"))
write_csv(spatial, paste0(home, "/output/spatial.csv"))
write_csv(spatial_density, paste0(home, "/output/spatial_density.csv"))
write_csv(range, paste0(home, "/output/range.csv"))
write_csv(genq, paste0(home, "/output/genq.csv"))
write_csv(index, paste0(home, "/output/index.csv"))
write_csv(dens_index, paste0(home, "/output/dens_index.csv"))

pars <- read_csv(paste0(home, "/output/pars.csv"))

Rows: 304 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): term, species
dbl (4): estimate, std.error, conf.low, conf.high

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

res <- read_csv(paste0(home, "/output/res.csv"))

Rows: 43773 Columns: 27
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (5): haul_id, species, period, life_stage, species_f
dbl (22): length_mm, number, year, day, month, lat, lon, temp_obs, yday, X, ...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

res_dens <- read_csv(paste0(home, "/output/res_dens.csv"))

Rows: 16456 Columns: 21
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): haul_id, species, period, national_haul_id
dbl (17): lat, lon, year, month, day, no_m2, yday, X, Y, temp, chl, depth, y...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

spatial <- read_csv(paste0(home, "/output/spatial.csv"))

Rows: 1270719 Columns: 20
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): species
dbl (19): X, Y, year, lon, lat, depth, temp, chl, temp_sc, temp_sq, chl_sc, ...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

spatial_density <- read_csv(paste0(home, "/output/spatial_density.csv"))

Rows: 1462968 Columns: 20
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): species
dbl (19): X, Y, year, lon, lat, depth, temp, chl, temp_sc, temp_sq, chl_sc, ...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

range <- read_csv(paste0(home, "/output/range.csv"))

Rows: 22 Columns: 8
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): term, model, species
dbl (5): estimate, std.error, conf.low, conf.high, cutoff

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

genq <- read_csv(paste0(home, "/output/genq.csv"))

Rows: 22 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): term, species
dbl (4): estimate, std.error, conf.low, conf.high

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

index <- read_csv(paste0(home, "/output/index.csv"))

Rows: 542 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): type, species
dbl (7): year, est, lwr, upr, log_est, se, pred_dens_sum

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

dens_index <- read_csv(paste0(home, "/output/dens_index.csv"))

Rows: 312 Columns: 8
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): type, species
dbl (6): year, est, lwr, upr, log_est, se

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Plot residuals and other things

# Plot residuals
ggplot(res, aes(observed, expected)) +
  geom_point(color = "grey30", shape = 21, size = 0.5) +
  geom_abline(col = "tomato3", linewidth = 0.6) +
  theme(aspect.ratio = 1) +
  labs(x = "Observed", y = "Expected") +
  facet_wrap(~species, scales = "free")

ggsave(paste0(home, "/figures/supp/residuals.pdf"), width = 20, height = 20, units = "cm")

# Density residuals
ggplot(res_dens, aes(observed, expected)) +
  geom_point(color = "grey30", shape = 21, size = 0.5) +
  geom_abline(col = "tomato3", linewidth = 0.6) +
  theme(aspect.ratio = 1) +
  labs(x = "Observed", y = "Expected") +
  facet_wrap(~species, scales = "free")

ggsave(paste0(home, "/figures/supp/residuals_density.pdf"), width = 20, height = 20, units = "cm")

# Plot range
ggplot(range, aes(estimate, reorder(species, estimate))) +
  geom_point() +
  geom_vline(xintercept = 5, alpha = 0.5, linetype = 2) +
  facet_wrap(~cutoff)

ggplot(range, aes(estimate / cutoff, reorder(species, estimate))) +
  geom_point() +
  geom_vline(xintercept = 2, alpha = 0.5, linetype = 2)

# Plot generalized gamma
# The lognormal occurs as the internal parameter gengamma_Q approaches 0
genq |>
  filter(term == "gengamma_Q") |>
  ggplot(aes(estimate, species)) +
  facet_wrap(~term, ncol = 1) +
  geom_point() +
  geom_vline(xintercept = 0, alpha = 0.5, linetype = 2)

filter: removed 11 rows (50%), 11 rows remaining

# If Q matches phi the distribution should be the gamma.
genq |>
  ggplot(aes(estimate, species, color = term)) +
  geom_point()

Coefficients

pars_order <- pars |>
  filter(term %in% c("temp_sc")) |>
  arrange(estimate)

pars |>
  # left_join(d |> distinct(life_stage, species)) |>
  # filter(life_stage == "larvae") |>
  filter(term %in% c("temp_sc", "chl_sc", "yday_ct")) |>
  mutate(life_stage = ifelse(species %in% c("Ammodytidae", "Syngnathus rostellatus", "Pomatoschistus sp.", 
                                            "Crystallogobius linearis", "Aphia minuta"), 
                             "Adult", "Larvae")) |> 
  mutate(term = fct_recode(term,
    "Day of the year" = "yday_ct",
    "Temperature" = "temp_sc",
    "Chlorophyll" = "chl_sc"
  )) |>
  mutate(
    sign = ifelse(estimate > 0, "pos", "neg"),
    sig = ifelse(estimate > 0 & conf.low > 0, "sig", "not sig"),
    sig = ifelse(estimate < 0 & conf.high < 0, "sig", sig)
  ) |>
  ggplot(aes(estimate, factor(species, levels = pars_order$species), fill = sig, color = sign, shape = sig)) +
  #facet_wrap(~ factor(term, levels = c("Temperature", "Day of the year", "Chlorophyll"))) +
  ggh4x::facet_grid2(life_stage ~ factor(term, levels = c("Temperature", "Day of the year", "Chlorophyll")),
                     scales = "free_y") +
  geom_point(fill = NA) +
  scale_shape_manual(values = c(21, 19)) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high),
    width = 0, alpha = 0.3
  ) +
  geom_vline(xintercept = 0, alpha = 0.3, linetype = 2, linewidth = 0.25) +
  theme(axis.text.y = element_text(face = "italic")) +
  labs(y = "Species", x = "Slope") +
  scale_color_brewer(palette = "Dark2", direction = -1) +
  scale_fill_brewer(palette = "Dark2", direction = -1) +
  guides(color = "none", shape = "none") +
  geom_stripped_rows(aes(y = species), inherit.aes = FALSE) +
  theme(legend.position = "bottom")

ggsave(paste0(home, "/figures/coefs.pdf"), width = 17, height = 10, units = "cm")

Spatial plots and indices

trends <- index |>
  filter(!species == "Sardina pilchardus") |> 
  filter(type == "weighted") |> 
  rename(ind = est, 
         ind_lwr = lwr,
         ind_upr = upr,
         ind_se = se) |> 
  group_by(species) |> 
  nest() |> 
  mutate(
    # Fit models
    model = map(data, ~sdmTMB(
      ind ~ s(year, k = 8),
      family = lognormal(),
      data = .x,
      spatial = "off"
    )),
    # Get predictions with CIs
    preds = map(model, ~predict(.x, se_fit = TRUE) |> 
      mutate(
        lwr = exp(est - 1.96 * est_se),
        upr = exp(est + 1.96 * est_se)
      ))
  ) |> 
  select(species, preds) |> 
  unnest(preds)

filter: no rows removed

filter: removed 271 rows (50%), 271 rows remaining

rename: renamed 4 variables (ind, ind_lwr, ind_upr, ind_se)

group_by: one grouping variable (species)

Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.
mutate: new variable 'lwr' (double) with 14 unique values and 0% NA

        new variable 'upr' (double) with 14 unique values and 0% NA

Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.
mutate: new variable 'lwr' (double) with 28 unique values and 0% NA

        new variable 'upr' (double) with 28 unique values and 0% NA

Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.
mutate: new variable 'lwr' (double) with 29 unique values and 0% NA

        new variable 'upr' (double) with 29 unique values and 0% NA

Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.
mutate: new variable 'lwr' (double) with 28 unique values and 0% NA

        new variable 'upr' (double) with 28 unique values and 0% NA

Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.
mutate: new variable 'lwr' (double) with 28 unique values and 0% NA

        new variable 'upr' (double) with 28 unique values and 0% NA

Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.
mutate: new variable 'lwr' (double) with 29 unique values and 0% NA

        new variable 'upr' (double) with 29 unique values and 0% NA

Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.
mutate: new variable 'lwr' (double) with 22 unique values and 0% NA

        new variable 'upr' (double) with 22 unique values and 0% NA

Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.
mutate: new variable 'lwr' (double) with 17 unique values and 0% NA

        new variable 'upr' (double) with 17 unique values and 0% NA

Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.
mutate: new variable 'lwr' (double) with 28 unique values and 0% NA

        new variable 'upr' (double) with 28 unique values and 0% NA

Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.
mutate: new variable 'lwr' (double) with 26 unique values and 0% NA

        new variable 'upr' (double) with 26 unique values and 0% NA

Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.
mutate: new variable 'lwr' (double) with 22 unique values and 0% NA

        new variable 'upr' (double) with 22 unique values and 0% NA

mutate (grouped): new variable 'model' (list) with 11 unique values and 0% NA

                  new variable 'preds' (list) with 11 unique values and 0% NA

select: dropped 2 variables (data, model)

ggplot(trends, aes(year, exp(est))) +
  facet_wrap(~species, scales = "free") +
  geom_point(aes(year, ind), alpha = 0.8, color = "steelblue") +
  geom_errorbar(aes(ymin = ind_lwr, ymax = ind_upr), color = "steelblue", width = 0, alpha = 0.5) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1) +
  geom_line(color = "tomato", alpha = 0.8) +
  theme(strip.text = element_text(face = "italic")) +
  labs(x = "Year", y = "Length (mm)")

ggsave(paste0(home, "/figures/size_index.pdf"), width = 21, height = 16, units = "cm")

trends |>
  filter(type == "weighted") |> 
  ungroup() |>
  mutate(year2 = ifelse(year == min(year), "min", NA), .by = species) |>
  mutate(year2 = ifelse(year == max(year), "max", year2), .by = species) |>
  filter(year2 %in% c("min", "max")) |>
  dplyr::select(year2, species, est) |>
  pivot_wider(values_from = est, names_from = year2) |>
  mutate(percent = ((`max` - `min`) / abs(`min`)) * 100) |>
  arrange(desc(percent)) |>
  ggplot(aes(percent, reorder(species, percent))) +
  geom_vline(xintercept = 0, alpha = 0.5, linetype = 2) +
  geom_point(color = "steelblue") +
  labs(x = "Percent change", y = "Species") +
  theme(axis.text.y = element_text(face = "italic"))

filter (grouped): no rows removed
ungroup: no grouping variables
mutate: new variable 'year2' (character) with 2 unique values and 96% NA
mutate: changed 11 values (4%) of 'year2' (11 fewer NA)
filter: removed 249 rows (92%), 22 rows remaining
pivot_wider: reorganized (year2, est) into (min, max) [was 22x3, now 11x3]
mutate: new variable 'percent' (double) with 11 unique values and 0% NA

# ggsave(paste0(home, "/figures/percent_change.pdf"), width = 12, height = 9, units = "cm")

# density index
trends_dens <- dens_index |>
  rename(ind = est, 
         ind_lwr = lwr,
         ind_upr = upr,
         ind_se = se) |> 
  group_by(species) |> 
  nest() |> 
  mutate(
    # Fit models
    model = map(data, ~sdmTMB(
      ind ~ s(year, k = 8),
      family = lognormal(),
      data = .x,
      spatial = "off"
    )),
    # Get predictions with CIs
    preds = map(model, ~predict(.x, se_fit = TRUE) |> 
      mutate(
        lwr = exp(est - 1.96 * est_se),
        upr = exp(est + 1.96 * est_se)
      ))
  ) |> 
  select(species, preds) |> 
  unnest(preds)

rename: renamed 4 variables (ind, ind_lwr, ind_upr, ind_se)
group_by: one grouping variable (species)
Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.mutate: new variable 'lwr' (double) with 27 unique values and 0% NA
        new variable 'upr' (double) with 27 unique values and 0% NA
Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.mutate: new variable 'lwr' (double) with 29 unique values and 0% NA
        new variable 'upr' (double) with 29 unique values and 0% NA
Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.mutate: new variable 'lwr' (double) with 29 unique values and 0% NA
        new variable 'upr' (double) with 29 unique values and 0% NA
Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.mutate: new variable 'lwr' (double) with 29 unique values and 0% NA
        new variable 'upr' (double) with 29 unique values and 0% NA
Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.mutate: new variable 'lwr' (double) with 29 unique values and 0% NA
        new variable 'upr' (double) with 29 unique values and 0% NA
Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.mutate: new variable 'lwr' (double) with 29 unique values and 0% NA
        new variable 'upr' (double) with 29 unique values and 0% NA
Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.mutate: new variable 'lwr' (double) with 28 unique values and 0% NA
        new variable 'upr' (double) with 28 unique values and 0% NA
Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.mutate: new variable 'lwr' (double) with 25 unique values and 0% NA
        new variable 'upr' (double) with 25 unique values and 0% NA
Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.mutate: new variable 'lwr' (double) with 29 unique values and 0% NA
        new variable 'upr' (double) with 29 unique values and 0% NA
Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.mutate: new variable 'lwr' (double) with 29 unique values and 0% NA
        new variable 'upr' (double) with 29 unique values and 0% NA
Prediction can be slow when `se_fit = TRUE` and random fields are included
(i.e., `re_form = NA`). Consider using the `nsim` argument to take draws from
the joint precision matrix and summarizing the standard devation of those
draws.mutate: new variable 'lwr' (double) with 29 unique values and 0% NA
        new variable 'upr' (double) with 29 unique values and 0% NA
mutate (grouped): new variable 'model' (list) with 11 unique values and 0% NA
                  new variable 'preds' (list) with 11 unique values and 0% NA
select: dropped 2 variables (data, model)

ggplot(trends_dens, aes(year, exp(est))) +
  facet_wrap(~species, scales = "free") +
  geom_point(aes(year, ind), alpha = 0.8, color = "steelblue") +
  geom_errorbar(aes(ymin = ind_lwr, ymax = ind_upr), color = "steelblue", width = 0, alpha = 0.5) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.1) +
  geom_line(color = "tomato", alpha = 0.8) +
  theme(strip.text = element_text(face = "italic")) +
  labs(x = "Year", y = "Relative abundance")

ggsave(paste0(home, "/figures/density_index.pdf"), width = 21, height = 16, units = "cm")

# Length
plot_map_fc +
  geom_raster(
    data = spatial |>
      summarise(est = mean(exp(est)), .by = c(X, Y, species)) |> 
      mutate(est_sc = est / mean(est), .by = species),
    aes(X * 1000, Y * 1000, fill = est_sc)
  ) +
  scale_fill_gradient2(midpoint = 1) +
  facet_wrap(~species) +
  labs(fill = "Scaled size") +
  theme(strip.text = element_text(face = "italic")) + 
  geom_sf()

summarise: now 51,579 rows and 4 columns, ungrouped
mutate: new variable 'est_sc' (double) with 51,579 unique values and 0% NA

ggsave(paste0(home, "/figures/spatial_size.pdf"), width = 17, height = 17, units = "cm")

# Density
plot_map_fc +
  geom_raster(
    data = spatial_density |>
      summarise(pred_dens = mean(pred_dens), .by = c(X, Y, species)) |> 
      mutate(est_sc = pred_dens / mean(pred_dens), .by = species),
    aes(X * 1000, Y * 1000, fill = est_sc)
  ) +
  #scale_fill_viridis(trans = "sqrt") +
  scale_fill_gradient2(midpoint = 1, trans = "sqrt") +
  facet_wrap(~species) +
  labs(fill = "Scaled density") +
  theme(strip.text = element_text(face = "italic")) + 
  geom_sf()

summarise: now 51,579 rows and 4 columns, ungrouped
mutate: new variable 'est_sc' (double) with 51,579 unique values and 0% NA

ggsave(paste0(home, "/figures/spatial_density.pdf"), width = 17, height = 17, units = "cm")

Percent change in index

trends_change <- trends |>
  filter(type == "weighted") |> 
  filter(year >= 2012) |> 
  ungroup() |>
  mutate(year2 = ifelse(year == min(year), "min", NA), .by = species) |>
  mutate(year2 = ifelse(year == max(year), "max", year2), .by = species) |>
  filter(year2 %in% c("min", "max")) |>
  mutate(est = exp(est)) |> 
  dplyr::select(year2, species, est) |>
  pivot_wider(values_from = est, names_from = year2) |>
  mutate(percent = ((`max` - `min`) / abs(`min`)) * 100) |>
  arrange(desc(percent))

filter (grouped): no rows removed

filter (grouped): removed 165 rows (61%), 106 rows remaining

ungroup: no grouping variables

mutate: new variable 'year2' (character) with 2 unique values and 90% NA

mutate: changed 11 values (10%) of 'year2' (11 fewer NA)

filter: removed 84 rows (79%), 22 rows remaining

mutate: changed 22 values (100%) of 'est' (0 new NA)

pivot_wider: reorganized (year2, est) into (min, max) [was 22x3, now 11x3]

mutate: new variable 'percent' (double) with 11 unique values and 0% NA

trends_change

# A tibble: 11 × 4
   species                    min   max percent
   <chr>                    <dbl> <dbl>   <dbl>
 1 Pholis gunnellus         13.9  14.2    2.03 
 2 Agonus cataphractus       8.97  8.99   0.212
 3 Aphia minuta             33.4  32.7   -2.11 
 4 Chirolophis ascanii      15.6  15.0   -4.37 
 5 Ammodytidae              32.5  31.0   -4.68 
 6 Clupea harengus          25.9  24.5   -5.52 
 7 Pomatoschistus sp        42.2  39.3   -6.98 
 8 Crystallogobius linearis 21.3  19.8   -7.32 
 9 Myoxocephalus scorpius    8.17  7.35  -9.98 
10 Syngnathus rostellatus   70.1  55.0  -21.5  
11 Microstomus kitt         24.3  18.4  -24.5  

trends_change |> 
  summarise(mean = mean(percent),
            max = max(percent),
            min = min(percent),
            median = median(percent))

summarise: now one row and 4 columns, ungrouped

# A tibble: 1 × 4
   mean   max   min median
  <dbl> <dbl> <dbl>  <dbl>
1 -7.70  2.03 -24.5  -5.52

trends_dens_change <- trends_dens |>
  filter(year >= 2012) |> 
  ungroup() |>
  mutate(year2 = ifelse(year == min(year), "min", NA), .by = species) |>
  mutate(year2 = ifelse(year == max(year), "max", year2), .by = species) |>
  filter(year2 %in% c("min", "max")) |>
  mutate(est = exp(est)) |> 
  dplyr::select(year2, species, est) |>
  pivot_wider(values_from = est, names_from = year2) |>
  mutate(percent = ((`max` - `min`) / abs(`min`)) * 100) |>
  arrange(desc(percent)) 

filter (grouped): removed 192 rows (62%), 120 rows remaining

ungroup: no grouping variables

mutate: new variable 'year2' (character) with 2 unique values and 91% NA

mutate: changed 11 values (9%) of 'year2' (11 fewer NA)

filter: removed 98 rows (82%), 22 rows remaining

mutate: changed 22 values (100%) of 'est' (0 new NA)

pivot_wider: reorganized (year2, est) into (min, max) [was 22x3, now 11x3]

mutate: new variable 'percent' (double) with 11 unique values and 0% NA

trends_dens_change

# A tibble: 11 × 4
   species                          min         max percent
   <chr>                          <dbl>       <dbl>   <dbl>
 1 Myoxocephalus scorpius     97071127.  118325921.   21.9 
 2 Ammodytidae               236250333.  285803937.   21.0 
 3 Microstomus kitt          153737814.  161623726.    5.13
 4 Chirolophis ascanii       162758888.  167447452.    2.88
 5 Pholis gunnellus          680394400.  597100062.  -12.2 
 6 Clupea harengus          1210993076.  999631894.  -17.5 
 7 Crystallogobius linearis 3626289793. 2897647849.  -20.1 
 8 Agonus cataphractus        63946975.   49074732.  -23.3 
 9 Syngnathus rostellatus    312917153.  224403920.  -28.3 
10 Aphia minuta             1217931138.  513127608.  -57.9 
11 Pomatoschistus sp         196607581.   78458813.  -60.1 

trends_dens_change |> 
  summarise(mean = mean(percent),
            max = max(percent),
            min = min(percent),
            median = median(percent))

summarise: now one row and 4 columns, ungrouped

# A tibble: 1 × 4
   mean   max   min median
  <dbl> <dbl> <dbl>  <dbl>
1 -15.3  21.9 -60.1  -17.5